Nektar++
VCSMapping.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File VCSMapping.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Velocity Correction Scheme w/ coordinate transformation
32 // for the Incompressible Navier Stokes equations
33 ///////////////////////////////////////////////////////////////////////////////
34 
37 #include <SolverUtils/Core/Misc.h>
38 
39 #include <boost/algorithm/string.hpp>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 string VCSMapping::className =
47  "VCSMapping", VCSMapping::create);
48 
49 /**
50  * Constructor. Creates ...
51  *
52  * \param
53  * \param
54  */
55 VCSMapping::VCSMapping(const LibUtilities::SessionReaderSharedPtr &pSession,
57  : UnsteadySystem(pSession, pGraph),
58  VelocityCorrectionScheme(pSession, pGraph)
59 {
60 }
61 
62 void VCSMapping::v_InitObject(bool DeclareField)
63 {
65 
67  ASSERTL0(m_mapping, "Could not create mapping in VCSMapping.");
68 
69  std::string vExtrapolation = "Mapping";
71  vExtrapolation, m_session, m_fields, m_pressure, m_velocity,
72  m_advObject);
73  m_extrapolation->SubSteppingTimeIntegration(m_intScheme);
74  m_extrapolation->GenerateHOPBCMap(m_session);
75 
76  // Storage to extrapolate pressure forcing
77  int physTot = m_fields[0]->GetTotPoints();
78  int intSteps = 1;
79 
80  if (m_intScheme->GetName() == "IMEX" ||
81  m_intScheme->GetName() == "IMEXGear")
82  {
83  m_intSteps = m_intScheme->GetOrder();
84  }
85  else
86  {
88  "Integration method not suitable: "
89  "Options include IMEXGear or IMEXOrder{1,2,3,4}");
90  }
91 
93  for (int i = 0; i < m_presForcingCorrection.size(); i++)
94  {
96  }
97  m_verbose = (m_session->DefinesCmdLineArgument("verbose")) ? true : false;
98 
99  // Load solve parameters related to the mapping
100  // Flags determining if pressure/viscous terms should be treated implicitly
101  m_session->MatchSolverInfo("MappingImplicitPressure", "True",
102  m_implicitPressure, false);
103  m_session->MatchSolverInfo("MappingImplicitViscous", "True",
104  m_implicitViscous, false);
105  m_session->MatchSolverInfo("MappingNeglectViscous", "True",
106  m_neglectViscous, false);
107 
108  if (m_neglectViscous)
109  {
110  m_implicitViscous = false;
111  }
112 
113  // Tolerances and relaxation parameters for implicit terms
114  m_session->LoadParameter("MappingPressureTolerance", m_pressureTolerance,
115  1e-12);
116  m_session->LoadParameter("MappingViscousTolerance", m_viscousTolerance,
117  1e-12);
118  m_session->LoadParameter("MappingPressureRelaxation", m_pressureRelaxation,
119  1.0);
120  m_session->LoadParameter("MappingViscousRelaxation", m_viscousRelaxation,
121  1.0);
122 }
123 
124 /**
125  * Destructor
126  */
128 {
129 }
130 
132 {
134 
135  // Set up Field Meta Data for output files
136  m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
137  m_fieldMetaDataMap["TimeStep"] =
138  boost::lexical_cast<std::string>(m_timestep);
139 
140  // Correct Dirichlet boundary conditions to account for mapping
141  m_mapping->UpdateBCs(0.0);
142  //
144  for (int i = 0; i < m_nConvectiveFields; ++i)
145  {
146  m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
147  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
148  m_fields[i]->UpdatePhys());
150  }
151 
152  // Initialise m_gradP
153  int physTot = m_fields[0]->GetTotPoints();
155  for (int i = 0; i < m_nConvectiveFields; ++i)
156  {
157  m_gradP[i] = Array<OneD, NekDouble>(physTot, 0.0);
159  m_pressure->GetPhys(), m_gradP[i]);
160  if (m_pressure->GetWaveSpace())
161  {
162  m_pressure->HomogeneousBwdTrans(m_gradP[i], m_gradP[i]);
163  }
164  }
165 }
166 
167 /**
168  * Explicit part of the method - Advection, Forcing + HOPBCs
169  */
171  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
172  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
173 {
174  EvaluateAdvectionTerms(inarray, outarray, time);
175 
176  // Smooth advection
177  if (m_SmoothAdvection)
178  {
179  for (int i = 0; i < m_nConvectiveFields; ++i)
180  {
181  m_pressure->SmoothField(outarray[i]);
182  }
183  }
184 
185  // Add forcing terms
186  for (auto &x : m_forcing)
187  {
188  x->Apply(m_fields, inarray, outarray, time);
189  }
190 
191  // Add mapping terms
192  ApplyIncNSMappingForcing(inarray, outarray);
193 
194  // Calculate High-Order pressure boundary conditions
195  m_extrapolation->EvaluatePressureBCs(inarray, outarray, m_kinvis);
196 
197  // Update mapping and deal with Dirichlet boundary conditions
198  if (m_mapping->IsTimeDependent())
199  {
200  if (m_mapping->IsFromFunction())
201  {
202  // If the transformation is explicitly defined, update it here
203  // Otherwise, it will be done somewhere else (ForcingMovingBody)
204  m_mapping->UpdateMapping(time + m_timestep);
205  }
206  m_mapping->UpdateBCs(time + m_timestep);
207  }
208 }
209 
210 /**
211  * Forcing term for Poisson solver solver
212  */
214  const Array<OneD, const Array<OneD, NekDouble>> &fields,
216 {
217  if (m_mapping->HasConstantJacobian())
218  {
220  aii_Dt);
221  }
222  else
223  {
224  int physTot = m_fields[0]->GetTotPoints();
225  int nvel = m_nConvectiveFields;
226  Array<OneD, NekDouble> wk(physTot, 0.0);
227 
228  Array<OneD, NekDouble> Jac(physTot, 0.0);
229  m_mapping->GetJacobian(Jac);
230 
231  // Calculate div(J*u/Dt)
232  Vmath::Zero(physTot, Forcing[0], 1);
233  for (int i = 0; i < nvel; ++i)
234  {
235  if (m_fields[i]->GetWaveSpace())
236  {
237  m_fields[i]->HomogeneousBwdTrans(fields[i], wk);
238  }
239  else
240  {
241  Vmath::Vcopy(physTot, fields[i], 1, wk, 1);
242  }
243  Vmath::Vmul(physTot, wk, 1, Jac, 1, wk, 1);
244  if (m_fields[i]->GetWaveSpace())
245  {
246  m_fields[i]->HomogeneousFwdTrans(wk, wk);
247  }
248  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[i], wk, wk);
249  Vmath::Vadd(physTot, wk, 1, Forcing[0], 1, Forcing[0], 1);
250  }
251  Vmath::Smul(physTot, 1.0 / aii_Dt, Forcing[0], 1, Forcing[0], 1);
252 
253  //
254  // If the mapping viscous terms are being treated explicitly
255  // we need to apply a correction to the forcing
256  if (!m_implicitViscous)
257  {
258  bool wavespace = m_fields[0]->GetWaveSpace();
259  m_fields[0]->SetWaveSpace(false);
260 
261  //
262  // Part 1: div(J*grad(U/J . grad(J)))
264  Array<OneD, Array<OneD, NekDouble>> velocity(nvel);
265  for (int i = 0; i < tmp.size(); i++)
266  {
267  tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
268  velocity[i] = Array<OneD, NekDouble>(physTot, 0.0);
269  if (wavespace)
270  {
271  m_fields[0]->HomogeneousBwdTrans(m_fields[i]->GetPhys(),
272  velocity[i]);
273  }
274  else
275  {
276  Vmath::Vcopy(physTot, m_fields[i]->GetPhys(), 1,
277  velocity[i], 1);
278  }
279  }
280  // Calculate wk = U.grad(J)
281  m_mapping->DotGradJacobian(velocity, wk);
282  // Calculate wk = (U.grad(J))/J
283  Vmath::Vdiv(physTot, wk, 1, Jac, 1, wk, 1);
284  // J*grad[(U.grad(J))/J]
285  for (int i = 0; i < nvel; ++i)
286  {
287  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], wk,
288  tmp[i]);
289  Vmath::Vmul(physTot, Jac, 1, tmp[i], 1, tmp[i], 1);
290  }
291  // div(J*grad[(U.grad(J))/J])
292  Vmath::Zero(physTot, wk, 1);
293  for (int i = 0; i < nvel; ++i)
294  {
295  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], tmp[i],
296  tmp[i]);
297  Vmath::Vadd(physTot, wk, 1, tmp[i], 1, wk, 1);
298  }
299 
300  // Part 2: grad(J) . curl(curl(U))
301  m_mapping->CurlCurlField(velocity, tmp, m_implicitViscous);
302  // dont need velocity any more, so reuse it
303  m_mapping->DotGradJacobian(tmp, velocity[0]);
304 
305  // Add two parts
306  Vmath::Vadd(physTot, velocity[0], 1, wk, 1, wk, 1);
307 
308  // Multiply by kinvis and prepare to extrapolate
309  int nlevels = m_presForcingCorrection.size();
310  Vmath::Smul(physTot, m_kinvis, wk, 1,
311  m_presForcingCorrection[nlevels - 1], 1);
312 
313  // Extrapolate correction
314  m_extrapolation->ExtrapolateArray(m_presForcingCorrection);
315 
316  // Put in wavespace
317  if (wavespace)
318  {
319  m_fields[0]->HomogeneousFwdTrans(
320  m_presForcingCorrection[nlevels - 1], wk);
321  }
322  else
323  {
324  Vmath::Vcopy(physTot, m_presForcingCorrection[nlevels - 1], 1,
325  wk, 1);
326  }
327  // Apply correction: Forcing = Forcing - correction
328  Vmath::Vsub(physTot, Forcing[0], 1, wk, 1, Forcing[0], 1);
329 
330  m_fields[0]->SetWaveSpace(wavespace);
331  }
332  }
333 }
334 
335 /**
336  * Forcing term for Helmholtz solver
337  */
339  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
341 {
342  NekDouble aii_dtinv = 1.0 / aii_Dt;
343  int physTot = m_fields[0]->GetTotPoints();
344 
345  // Grad p
346  m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
347 
348  int nvel = m_velocity.size();
349  if (nvel == 2)
350  {
351  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1]);
352  }
353  else
354  {
355  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1],
356  Forcing[2]);
357  }
358 
359  // Copy grad p in physical space to m_gradP to reuse later
360  if (m_pressure->GetWaveSpace())
361  {
362  for (int i = 0; i < nvel; i++)
363  {
364  m_pressure->HomogeneousBwdTrans(Forcing[i], m_gradP[i]);
365  }
366  }
367  else
368  {
369  for (int i = 0; i < nvel; i++)
370  {
371  Vmath::Vcopy(physTot, Forcing[i], 1, m_gradP[i], 1);
372  }
373  }
374 
375  if ((!m_mapping->HasConstantJacobian()) || m_implicitPressure)
376  {
377  // If pressure terms are treated explicitly, we need to divide by J
378  // if they are implicit, we need to calculate G(p)
379  if (m_implicitPressure)
380  {
381  m_mapping->RaiseIndex(m_gradP, Forcing);
382  }
383  else
384  {
385  Array<OneD, NekDouble> Jac(physTot, 0.0);
386  m_mapping->GetJacobian(Jac);
387  for (int i = 0; i < nvel; i++)
388  {
389  Vmath::Vdiv(physTot, m_gradP[i], 1, Jac, 1, Forcing[i], 1);
390  }
391  }
392  // Transform back to wavespace
393  if (m_pressure->GetWaveSpace())
394  {
395  for (int i = 0; i < nvel; i++)
396  {
397  m_pressure->HomogeneousFwdTrans(Forcing[i], Forcing[i]);
398  }
399  }
400  }
401 
402  // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
403  // need to be updated for the convected fields.
404  for (int i = 0; i < nvel; ++i)
405  {
406  Blas::Daxpy(physTot, -aii_dtinv, inarray[i], 1, Forcing[i], 1);
407  Blas::Dscal(physTot, 1.0 / m_kinvis, &(Forcing[i])[0], 1);
408  }
409 }
410 
411 /**
412  * Solve pressure system
413  */
415 {
416  if (!m_implicitPressure)
417  {
419  }
420  else
421  {
422  int physTot = m_fields[0]->GetTotPoints();
423  int nvel = m_nConvectiveFields;
424  bool converged = false; // flag to mark if system converged
425  int s = 0; // iteration counter
426  NekDouble error; // L2 error at current iteration
427  NekDouble forcing_L2 = 0.0; // L2 norm of F
428 
429  int maxIter;
430  m_session->LoadParameter("MappingMaxIter", maxIter, 5000);
431 
432  // rhs of the equation at current iteration
433  Array<OneD, NekDouble> F_corrected(physTot, 0.0);
434  // Pressure field at previous iteration
435  Array<OneD, NekDouble> previous_iter(physTot, 0.0);
436  // Temporary variables
440  for (int i = 0; i < nvel; ++i)
441  {
442  wk1[i] = Array<OneD, NekDouble>(physTot, 0.0);
443  wk2[i] = Array<OneD, NekDouble>(physTot, 0.0);
444  gradP[i] = Array<OneD, NekDouble>(physTot, 0.0);
445  }
446 
447  // Jacobian
448  Array<OneD, NekDouble> Jac(physTot, 0.0);
449  m_mapping->GetJacobian(Jac);
450 
451  // Factors for Laplacian system
453  factors[StdRegions::eFactorLambda] = 0.0;
454 
455  m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
456  forcing_L2 = m_pressure->L2(Forcing, wk1[0]);
457  while (!converged)
458  {
459  // Update iteration counter and set previous iteration field
460  // (use previous timestep solution for first iteration)
461  s++;
462  ASSERTL0(s < maxIter,
463  "VCSMapping exceeded maximum number of iterations.");
464 
465  Vmath::Vcopy(physTot, m_pressure->GetPhys(), 1, previous_iter, 1);
466 
467  // Correct pressure bc to account for iteration
468  m_extrapolation->CorrectPressureBCs(previous_iter);
469 
470  //
471  // Calculate forcing term for this iteration
472  //
473  for (int i = 0; i < nvel; ++i)
474  {
476  previous_iter, gradP[i]);
477  if (m_pressure->GetWaveSpace())
478  {
479  m_pressure->HomogeneousBwdTrans(gradP[i], wk1[i]);
480  }
481  else
482  {
483  Vmath::Vcopy(physTot, gradP[i], 1, wk1[i], 1);
484  }
485  }
486  m_mapping->RaiseIndex(wk1, wk2); // G(p)
487 
488  m_mapping->Divergence(wk2, F_corrected); // div(G(p))
489  if (!m_mapping->HasConstantJacobian())
490  {
491  Vmath::Vmul(physTot, F_corrected, 1, Jac, 1, F_corrected, 1);
492  }
493  // alpha*J*div(G(p))
494  Vmath::Smul(physTot, m_pressureRelaxation, F_corrected, 1,
495  F_corrected, 1);
496  if (m_pressure->GetWaveSpace())
497  {
498  m_pressure->HomogeneousFwdTrans(F_corrected, F_corrected);
499  }
500  // alpha*J*div(G(p)) - p_ii
501  for (int i = 0; i < m_nConvectiveFields; ++i)
502  {
504  gradP[i], wk1[0]);
505  Vmath::Vsub(physTot, F_corrected, 1, wk1[0], 1, F_corrected, 1);
506  }
507  // p_i,i - J*div(G(p))
508  Vmath::Neg(physTot, F_corrected, 1);
509  // alpha*F - alpha*J*div(G(p)) + p_i,i
510  Vmath::Smul(physTot, m_pressureRelaxation, Forcing, 1, wk1[0], 1);
511  Vmath::Vadd(physTot, wk1[0], 1, F_corrected, 1, F_corrected, 1);
512 
513  //
514  // Solve system
515  //
516  m_pressure->HelmSolve(F_corrected, m_pressure->UpdateCoeffs(),
517  factors);
518  m_pressure->BwdTrans(m_pressure->GetCoeffs(),
519  m_pressure->UpdatePhys());
520 
521  //
522  // Test convergence
523  //
524  error = m_pressure->L2(m_pressure->GetPhys(), previous_iter);
525  if (forcing_L2 != 0)
526  {
527  if ((error / forcing_L2 < m_pressureTolerance))
528  {
529  converged = true;
530  }
531  }
532  else
533  {
534  if (error < m_pressureTolerance)
535  {
536  converged = true;
537  }
538  }
539  }
540  if (m_verbose && m_session->GetComm()->GetRank() == 0)
541  {
542  std::cout << " Pressure system (mapping) converged in " << s
543  << " iterations with error = " << error << std::endl;
544  }
545  }
546 }
547 
548 /**
549  * Solve velocity system
550  */
552  const Array<OneD, const Array<OneD, NekDouble>> &Forcing,
553  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
554  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble aii_Dt)
555 {
556  boost::ignore_unused(inarray);
557 
558  if (!m_implicitViscous)
559  {
561  aii_Dt);
562  }
563  else
564  {
565  int physTot = m_fields[0]->GetTotPoints();
566  int nvel = m_nConvectiveFields;
567  bool converged = false; // flag to mark if system converged
568  int s = 0; // iteration counter
569  NekDouble error, max_error; // L2 error at current iteration
570 
571  int maxIter;
572  m_session->LoadParameter("MappingMaxIter", maxIter, 5000);
573 
574  // L2 norm of F
576 
577  // rhs of the equation at current iteration
578  Array<OneD, Array<OneD, NekDouble>> F_corrected(nvel);
579  // Solution at previous iteration
580  Array<OneD, Array<OneD, NekDouble>> previous_iter(nvel);
581  // Working space
583  for (int i = 0; i < nvel; ++i)
584  {
585  F_corrected[i] = Array<OneD, NekDouble>(physTot, 0.0);
586  previous_iter[i] = Array<OneD, NekDouble>(physTot, 0.0);
587  wk[i] = Array<OneD, NekDouble>(physTot, 0.0);
588  }
589 
590  // Factors for Helmholtz system
592  factors[StdRegions::eFactorLambda] =
593  1.0 * m_viscousRelaxation / aii_Dt / m_kinvis;
594  if (m_useSpecVanVisc)
595  {
599  }
600 
601  // Calculate L2-norm of F and set initial solution for iteration
602  for (int i = 0; i < nvel; ++i)
603  {
604  forcing_L2[i] = m_fields[0]->L2(Forcing[i], wk[0]);
605  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), previous_iter[i]);
606  }
607 
608  while (!converged)
609  {
610  converged = true;
611  // Iteration counter
612  s++;
613  ASSERTL0(s < maxIter,
614  "VCSMapping exceeded maximum number of iterations.");
615 
616  max_error = 0.0;
617 
618  //
619  // Calculate forcing term for next iteration
620  //
621 
622  // Calculate L(U)- in this parts all components might be coupled
623  if (m_fields[0]->GetWaveSpace())
624  {
625  for (int i = 0; i < nvel; ++i)
626  {
627  m_fields[0]->HomogeneousBwdTrans(previous_iter[i], wk[i]);
628  }
629  }
630  else
631  {
632  for (int i = 0; i < nvel; ++i)
633  {
634  Vmath::Vcopy(physTot, previous_iter[i], 1, wk[i], 1);
635  }
636  }
637 
638  // (L(U^i) - 1/alpha*U^i_jj)
639  m_mapping->VelocityLaplacian(wk, F_corrected,
640  1.0 / m_viscousRelaxation);
641 
642  if (m_fields[0]->GetWaveSpace())
643  {
644  for (int i = 0; i < nvel; ++i)
645  {
646  m_fields[0]->HomogeneousFwdTrans(F_corrected[i],
647  F_corrected[i]);
648  }
649  }
650  else
651  {
652  for (int i = 0; i < nvel; ++i)
653  {
654  Vmath::Vcopy(physTot, F_corrected[i], 1, F_corrected[i], 1);
655  }
656  }
657 
658  // Loop velocity components
659  for (int i = 0; i < nvel; ++i)
660  {
661  // (-alpha*L(U^i) + U^i_jj)
662  Vmath::Smul(physTot, -1.0 * m_viscousRelaxation, F_corrected[i],
663  1, F_corrected[i], 1);
664  // F_corrected = alpha*F + (-alpha*L(U^i) + U^i_jj)
665  Vmath::Smul(physTot, m_viscousRelaxation, Forcing[i], 1, wk[0],
666  1);
667  Vmath::Vadd(physTot, wk[0], 1, F_corrected[i], 1,
668  F_corrected[i], 1);
669 
670  //
671  // Solve System
672  //
673  m_fields[i]->HelmSolve(F_corrected[i],
674  m_fields[i]->UpdateCoeffs(), factors);
675  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
676 
677  //
678  // Test convergence
679  //
680  error = m_fields[i]->L2(outarray[i], previous_iter[i]);
681 
682  if (forcing_L2[i] != 0)
683  {
684  if ((error / forcing_L2[i] >= m_viscousTolerance))
685  {
686  converged = false;
687  }
688  }
689  else
690  {
691  if (error >= m_viscousTolerance)
692  {
693  converged = false;
694  }
695  }
696  if (error > max_error)
697  {
698  max_error = error;
699  }
700 
701  // Copy field to previous_iter
702  Vmath::Vcopy(physTot, outarray[i], 1, previous_iter[i], 1);
703  }
704  }
705  if (m_verbose && m_session->GetComm()->GetRank() == 0)
706  {
707  std::cout << " Velocity system (mapping) converged in " << s
708  << " iterations with error = " << max_error << std::endl;
709  }
710  }
711 }
712 
713 /**
714  * Explicit terms of the mapping
715  */
717  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
718  Array<OneD, Array<OneD, NekDouble>> &outarray)
719 {
720  int physTot = m_fields[0]->GetTotPoints();
725  for (int i = 0; i < m_nConvectiveFields; ++i)
726  {
727  velPhys[i] = Array<OneD, NekDouble>(physTot, 0.0);
728  Forcing[i] = Array<OneD, NekDouble>(physTot, 0.0);
729  tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
730  }
731 
732  // Get fields and store velocity in wavespace and physical space
733  if (m_fields[0]->GetWaveSpace())
734  {
735  for (int i = 0; i < m_nConvectiveFields; ++i)
736  {
737  vel[i] = inarray[i];
738  m_fields[0]->HomogeneousBwdTrans(vel[i], velPhys[i]);
739  }
740  }
741  else
742  {
743  for (int i = 0; i < m_nConvectiveFields; ++i)
744  {
745  vel[i] = inarray[i];
746  Vmath::Vcopy(physTot, inarray[i], 1, velPhys[i], 1);
747  }
748  }
749 
750  // Advection contribution
752 
753  // Time-derivative contribution
754  if (m_mapping->IsTimeDependent())
755  {
756  MappingAccelerationCorrection(vel, velPhys, tmp);
757  for (int i = 0; i < m_nConvectiveFields; ++i)
758  {
759  Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
760  }
761  }
762 
763  // Pressure contribution
764  if (!m_implicitPressure)
765  {
767  for (int i = 0; i < m_nConvectiveFields; ++i)
768  {
769  Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
770  }
771  }
772  // Viscous contribution
773  if ((!m_implicitViscous) && (!m_neglectViscous))
774  {
775  MappingViscousCorrection(velPhys, tmp);
776  for (int i = 0; i < m_nConvectiveFields; ++i)
777  {
778  Vmath::Smul(physTot, m_kinvis, tmp[i], 1, tmp[i], 1);
779  Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
780  }
781  }
782 
783  // If necessary, transform to wavespace
784  if (m_fields[0]->GetWaveSpace())
785  {
786  for (int i = 0; i < m_nConvectiveFields; ++i)
787  {
788  m_fields[0]->HomogeneousFwdTrans(Forcing[i], Forcing[i]);
789  }
790  }
791 
792  // Add to outarray
793  for (int i = 0; i < m_nConvectiveFields; ++i)
794  {
795  Vmath::Vadd(physTot, outarray[i], 1, Forcing[i], 1, outarray[i], 1);
796  }
797 }
798 
800  const Array<OneD, const Array<OneD, NekDouble>> &velPhys,
801  Array<OneD, Array<OneD, NekDouble>> &outarray)
802 {
803  int physTot = m_fields[0]->GetTotPoints();
804  int nvel = m_nConvectiveFields;
805 
806  Array<OneD, Array<OneD, NekDouble>> wk(nvel * nvel);
807 
808  // Apply Christoffel symbols to obtain {i,kj}vel(k)
809  m_mapping->ApplyChristoffelContravar(velPhys, wk);
810 
811  // Calculate correction -U^j*{i,kj}vel(k)
812  for (int i = 0; i < nvel; i++)
813  {
814  Vmath::Zero(physTot, outarray[i], 1);
815  for (int j = 0; j < nvel; j++)
816  {
817  Vmath::Vvtvp(physTot, wk[i * nvel + j], 1, velPhys[j], 1,
818  outarray[i], 1, outarray[i], 1);
819  }
820  Vmath::Neg(physTot, outarray[i], 1);
821  }
822 }
823 
825  const Array<OneD, const Array<OneD, NekDouble>> &vel,
826  const Array<OneD, const Array<OneD, NekDouble>> &velPhys,
827  Array<OneD, Array<OneD, NekDouble>> &outarray)
828 {
829  int physTot = m_fields[0]->GetTotPoints();
830  int nvel = m_nConvectiveFields;
831 
832  Array<OneD, Array<OneD, NekDouble>> wk(nvel * nvel);
834  Array<OneD, Array<OneD, NekDouble>> coordVel(nvel);
835  for (int i = 0; i < nvel; i++)
836  {
837  tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
838  coordVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
839  }
840  // Get coordinates velocity in transformed system
841  m_mapping->GetCoordVelocity(tmp);
842  m_mapping->ContravarFromCartesian(tmp, coordVel);
843 
844  // Calculate first term: U^j u^i,j = U^j (du^i/dx^j + {i,kj}u^k)
845  m_mapping->ApplyChristoffelContravar(velPhys, wk);
846  for (int i = 0; i < nvel; i++)
847  {
848  Vmath::Zero(physTot, outarray[i], 1);
849 
850  m_fields[0]->PhysDeriv(velPhys[i], tmp[0], tmp[1]);
851  for (int j = 0; j < nvel; j++)
852  {
853  if (j == 2)
854  {
855  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j], vel[i],
856  tmp[2]);
857  if (m_fields[0]->GetWaveSpace())
858  {
859  m_fields[0]->HomogeneousBwdTrans(tmp[2], tmp[2]);
860  }
861  }
862 
863  Vmath::Vadd(physTot, wk[i * nvel + j], 1, tmp[j], 1,
864  wk[i * nvel + j], 1);
865 
866  Vmath::Vvtvp(physTot, coordVel[j], 1, wk[i * nvel + j], 1,
867  outarray[i], 1, outarray[i], 1);
868  }
869  }
870 
871  // Set wavespace to false and store current value
872  bool wavespace = m_fields[0]->GetWaveSpace();
873  m_fields[0]->SetWaveSpace(false);
874 
875  // Add -u^j U^i,j
876  m_mapping->ApplyChristoffelContravar(coordVel, wk);
877  for (int i = 0; i < nvel; i++)
878  {
879  if (nvel == 2)
880  {
881  m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1]);
882  }
883  else
884  {
885  m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1], tmp[2]);
886  }
887 
888  for (int j = 0; j < nvel; j++)
889  {
890  Vmath::Vadd(physTot, wk[i * nvel + j], 1, tmp[j], 1,
891  wk[i * nvel + j], 1);
892  Vmath::Neg(physTot, wk[i * nvel + j], 1);
893 
894  Vmath::Vvtvp(physTot, velPhys[j], 1, wk[i * nvel + j], 1,
895  outarray[i], 1, outarray[i], 1);
896  }
897  }
898 
899  // Restore value of wavespace
900  m_fields[0]->SetWaveSpace(wavespace);
901 }
902 
904  Array<OneD, Array<OneD, NekDouble>> &outarray)
905 {
906  int physTot = m_fields[0]->GetTotPoints();
907  int nvel = m_nConvectiveFields;
908 
909  // Calculate g^(ij)p_(,j)
910  m_mapping->RaiseIndex(m_gradP, outarray);
911 
912  // Calculate correction = (nabla p)/J - g^(ij)p_,j
913  // (Jac is not required if it is constant)
914  if (!m_mapping->HasConstantJacobian())
915  {
916  Array<OneD, NekDouble> Jac(physTot, 0.0);
917  m_mapping->GetJacobian(Jac);
918  for (int i = 0; i < nvel; ++i)
919  {
920  Vmath::Vdiv(physTot, m_gradP[i], 1, Jac, 1, m_gradP[i], 1);
921  }
922  }
923  for (int i = 0; i < nvel; ++i)
924  {
925  Vmath::Vsub(physTot, m_gradP[i], 1, outarray[i], 1, outarray[i], 1);
926  }
927 }
928 
930  const Array<OneD, const Array<OneD, NekDouble>> &velPhys,
931  Array<OneD, Array<OneD, NekDouble>> &outarray)
932 {
933  // L(U) - 1.0*d^2(u^i)/dx^jdx^j
934  m_mapping->VelocityLaplacian(velPhys, outarray, 1.0);
935 }
936 
937 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
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:269
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
NekDouble m_kinvis
Kinematic viscosity.
bool m_SmoothAdvection
bool to identify if advection term smoothing is requested
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
int m_intSteps
Number of time integration steps AND Order of extrapolation for pressure boundary conditions.
int m_nConvectiveFields
Number of fields to be convected;.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
NekDouble m_timestep
Time step size.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SOLVER_UTILS_EXPORT int GetTotPoints()
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:73
Base class for unsteady solvers.
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Sets up initial conditions.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
virtual void v_SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
Definition: VCSMapping.cpp:338
virtual ~VCSMapping()
Definition: VCSMapping.cpp:127
void MappingAccelerationCorrection(const Array< OneD, const Array< OneD, NekDouble >> &vel, const Array< OneD, const Array< OneD, NekDouble >> &velPhys, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: VCSMapping.cpp:824
void MappingViscousCorrection(const Array< OneD, const Array< OneD, NekDouble >> &velPhys, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: VCSMapping.cpp:929
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
Definition: VCSMapping.cpp:414
void ApplyIncNSMappingForcing(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: VCSMapping.cpp:716
NekDouble m_pressureTolerance
Definition: VCSMapping.h:86
NekDouble m_viscousRelaxation
Definition: VCSMapping.h:89
NekDouble m_pressureRelaxation
Definition: VCSMapping.h:88
NekDouble m_viscousTolerance
Definition: VCSMapping.h:87
void MappingPressureCorrection(Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: VCSMapping.cpp:903
GlobalMapping::MappingSharedPtr m_mapping
Definition: VCSMapping.h:75
void MappingAdvectionCorrection(const Array< OneD, const Array< OneD, NekDouble >> &velPhys, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: VCSMapping.cpp:799
virtual void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
Definition: VCSMapping.cpp:62
Array< OneD, Array< OneD, NekDouble > > m_presForcingCorrection
Definition: VCSMapping.h:117
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble >> &fields, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
Definition: VCSMapping.cpp:213
virtual void v_DoInitialise(void)
Sets up initial conditions.
Definition: VCSMapping.cpp:131
virtual void v_EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Definition: VCSMapping.cpp:170
virtual void v_SolveViscous(const Array< OneD, const Array< OneD, NekDouble >> &Forcing, const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble aii_Dt)
Definition: VCSMapping.cpp:551
Array< OneD, Array< OneD, NekDouble > > m_gradP
Definition: VCSMapping.h:92
NekDouble m_sVVCutoffRatio
cutt off ratio from which to start decayhing modes
NekDouble m_sVVDiffCoeff
Diffusion coefficient of SVV modes.
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
Array< OneD, Array< OneD, NekDouble > > m_F
virtual void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble >> &fields, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
virtual void v_SolveViscous(const Array< OneD, const Array< OneD, NekDouble >> &Forcing, const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble aii_Dt)
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:168
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:154
std::shared_ptr< SessionReader > SessionReaderSharedPtr
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:89
EquationSystemFactory & GetEquationSystemFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:48
double NekDouble
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:209
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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:574
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:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
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:284
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
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:419