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