Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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(m_session);
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  //
161  for(int i = 0; i < m_nConvectiveFields; ++i)
162  {
163  m_fields[i]->LocalToGlobal();
164  m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
165  m_fields[i]->GlobalToLocal();
166  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
167  m_fields[i]->UpdatePhys());
168  m_F[i] = Array< OneD, NekDouble> (m_fields[0]->GetTotPoints(), 0.0);
169  }
170 
171  // Initialise m_gradP
172  int physTot = m_fields[0]->GetTotPoints();
174  for(int i = 0; i < m_nConvectiveFields; ++i)
175  {
176  m_gradP[i] = Array<OneD, NekDouble>(physTot,0.0);
178  m_pressure->GetPhys(),
179  m_gradP[i]);
180  if(m_pressure->GetWaveSpace())
181  {
182  m_pressure->HomogeneousBwdTrans(m_gradP[i],m_gradP[i]);
183  }
184  }
185  }
186 
187  /**
188  * Explicit part of the method - Advection, Forcing + HOPBCs
189  */
191  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
192  Array<OneD, Array<OneD, NekDouble> > &outarray,
193  const NekDouble time)
194  {
195  EvaluateAdvectionTerms(inarray, outarray);
196 
197  // Smooth advection
199  {
200  for(int i = 0; i < m_nConvectiveFields; ++i)
201  {
202  m_pressure->SmoothField(outarray[i]);
203  }
204  }
205 
206  // Add forcing terms
207  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
208  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
209  {
210  (*x)->Apply(m_fields, inarray, outarray, time);
211  }
212 
213  // Add mapping terms
214  ApplyIncNSMappingForcing( inarray, outarray);
215 
216  // Calculate High-Order pressure boundary conditions
217  m_extrapolation->EvaluatePressureBCs(inarray,outarray,m_kinvis);
218 
219  // Update mapping and deal with Dirichlet boundary conditions
220  if (m_mapping->IsTimeDependent())
221  {
222  if (m_mapping->IsFromFunction())
223  {
224  // If the transformation is explicitly defined, update it here
225  // Otherwise, it will be done somewhere else (ForcingMovingBody)
226  m_mapping->UpdateMapping(time+m_timestep);
227  }
228  m_mapping->UpdateBCs(time+m_timestep);
229  }
230  }
231 
232 
233  /**
234  * Forcing term for Poisson solver solver
235  */
237  const Array<OneD, const Array<OneD, NekDouble> > &fields,
239  const NekDouble aii_Dt)
240  {
241  if (m_mapping->HasConstantJacobian())
242  {
244  Forcing, aii_Dt);
245  }
246  else
247  {
248  int physTot = m_fields[0]->GetTotPoints();
249  int nvel = m_nConvectiveFields;
250  Array<OneD, NekDouble> wk(physTot, 0.0);
251 
252  Array<OneD, NekDouble> Jac(physTot,0.0);
253  m_mapping->GetJacobian(Jac);
254 
255  // Calculate div(J*u/Dt)
256  Vmath::Zero(physTot,Forcing[0],1);
257  for(int i = 0; i < nvel; ++i)
258  {
259  if (m_fields[i]->GetWaveSpace())
260  {
261  m_fields[i]->HomogeneousBwdTrans(fields[i],wk);
262  }
263  else
264  {
265  Vmath::Vcopy(physTot, fields[i], 1, wk, 1);
266  }
267  Vmath::Vmul(physTot,wk,1,Jac,1,wk,1);
268  if (m_fields[i]->GetWaveSpace())
269  {
270  m_fields[i]->HomogeneousFwdTrans(wk,wk);
271  }
272  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[i],wk, wk);
273  Vmath::Vadd(physTot,wk,1,Forcing[0],1,Forcing[0],1);
274  }
275  Vmath::Smul(physTot,1.0/aii_Dt,Forcing[0],1,Forcing[0],1);
276 
277  //
278  // If the mapping viscous terms are being treated explicitly
279  // we need to apply a correction to the forcing
280  if (!m_implicitViscous)
281  {
282  bool wavespace = m_fields[0]->GetWaveSpace();
283  m_fields[0]->SetWaveSpace(false);
284 
285  //
286  // Part 1: div(J*grad(U/J . grad(J)))
288  Array<OneD, Array<OneD, NekDouble> > velocity (nvel);
289  for(int i = 0; i < tmp.num_elements(); i++)
290  {
291  tmp[i] = Array<OneD, NekDouble>(physTot,0.0);
292  velocity[i] = Array<OneD, NekDouble>(physTot,0.0);
293  if (wavespace)
294  {
295  m_fields[0]->HomogeneousBwdTrans(m_fields[i]->GetPhys(),
296  velocity[i]);
297  }
298  else
299  {
300  Vmath::Vcopy(physTot, m_fields[i]->GetPhys(), 1,
301  velocity[i], 1);
302  }
303  }
304  // Calculate wk = U.grad(J)
305  m_mapping->DotGradJacobian(velocity, wk);
306  // Calculate wk = (U.grad(J))/J
307  Vmath::Vdiv(physTot, wk, 1, Jac, 1, wk, 1);
308  // J*grad[(U.grad(J))/J]
309  for(int i = 0; i < nvel; ++i)
310  {
311  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i],
312  wk, tmp[i]);
313  Vmath::Vmul(physTot, Jac, 1, tmp[i], 1, tmp[i], 1);
314  }
315  // div(J*grad[(U.grad(J))/J])
316  Vmath::Zero(physTot, wk, 1);
317  for(int i = 0; i < nvel; ++i)
318  {
319  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i],
320  tmp[i], tmp[i]);
321  Vmath::Vadd(physTot, wk, 1, tmp[i], 1, wk, 1);
322  }
323 
324  // Part 2: grad(J) . curl(curl(U))
325  m_mapping->CurlCurlField(velocity, tmp, m_implicitViscous);
326  // dont need velocity any more, so reuse it
327  m_mapping->DotGradJacobian(tmp, velocity[0]);
328 
329  // Add two parts
330  Vmath::Vadd(physTot, velocity[0], 1, wk, 1, wk, 1);
331 
332  // Multiply by kinvis and prepare to extrapolate
333  int nlevels = m_presForcingCorrection.num_elements();
334  Vmath::Smul(physTot, m_kinvis, wk, 1,
335  m_presForcingCorrection[nlevels-1], 1);
336 
337  // Extrapolate correction
338  m_extrapolation->ExtrapolateArray(m_presForcingCorrection);
339 
340  // Put in wavespace
341  if (wavespace)
342  {
343  m_fields[0]->HomogeneousFwdTrans(
344  m_presForcingCorrection[nlevels-1],wk);
345  }
346  else
347  {
348  Vmath::Vcopy(physTot, m_presForcingCorrection[nlevels-1], 1,
349  wk, 1);
350  }
351  // Apply correction: Forcing = Forcing - correction
352  Vmath::Vsub(physTot, Forcing[0], 1, wk, 1, Forcing[0], 1);
353 
354  m_fields[0]->SetWaveSpace(wavespace);
355  }
356  }
357  }
358 
359  /**
360  * Forcing term for Helmholtz solver
361  */
363  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
365  const NekDouble aii_Dt)
366  {
367  NekDouble aii_dtinv = 1.0/aii_Dt;
368  int physTot = m_fields[0]->GetTotPoints();
369 
370  // Grad p
371  m_pressure->BwdTrans(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
372 
373  int nvel = m_velocity.num_elements();
374  if(nvel == 2)
375  {
376  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1]);
377  }
378  else
379  {
380  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1],
381  Forcing[2]);
382  }
383 
384  // Copy grad p in physical space to m_gradP to reuse later
385  if (m_pressure->GetWaveSpace())
386  {
387  for (int i=0; i<nvel; i++)
388  {
389  m_pressure->HomogeneousBwdTrans(Forcing[i],m_gradP[i]);
390  }
391  }
392  else
393  {
394  for (int i=0; i<nvel; i++)
395  {
396  Vmath::Vcopy(physTot, Forcing[i], 1, m_gradP[i], 1);
397  }
398  }
399 
400  if ( (!m_mapping->HasConstantJacobian()) || m_implicitPressure)
401  {
402  // If pressure terms are treated explicitly, we need to divide by J
403  // if they are implicit, we need to calculate G(p)
404  if (m_implicitPressure)
405  {
406  m_mapping->RaiseIndex(m_gradP, Forcing);
407  }
408  else
409  {
410  Array<OneD, NekDouble> Jac(physTot,0.0);
411  m_mapping->GetJacobian(Jac);
412  for (int i=0; i<nvel; i++)
413  {
414  Vmath::Vdiv(physTot, m_gradP[i], 1, Jac, 1, Forcing[i], 1);
415  }
416  }
417  // Transform back to wavespace
418  if (m_pressure->GetWaveSpace())
419  {
420  for (int i=0; i<nvel; i++)
421  {
422  m_pressure->HomogeneousFwdTrans(Forcing[i],Forcing[i]);
423  }
424  }
425  }
426 
427  // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
428  // need to be updated for the convected fields.
429  for(int i = 0; i < nvel; ++i)
430  {
431  Blas::Daxpy(physTot,-aii_dtinv,inarray[i],1,Forcing[i],1);
432  Blas::Dscal(physTot,1.0/m_kinvis,&(Forcing[i])[0],1);
433  }
434  }
435 
436  /**
437  * Solve pressure system
438  */
441  {
442  if (!m_implicitPressure)
443  {
445  }
446  else
447  {
448  int physTot = m_fields[0]->GetTotPoints();
449  int nvel = m_nConvectiveFields;
450  bool converged = false; // flag to mark if system converged
451  int s = 0; // iteration counter
452  NekDouble error; // L2 error at current iteration
453  NekDouble forcing_L2 = 0.0; // L2 norm of F
454 
455  int maxIter;
456  m_session->LoadParameter("MappingMaxIter",maxIter,5000);
457 
458  // rhs of the equation at current iteration
459  Array< OneD, NekDouble> F_corrected(physTot, 0.0);
460  // Pressure field at previous iteration
461  Array<OneD, NekDouble> previous_iter (physTot, 0.0);
462  // Temporary variables
466  for(int i = 0; i < nvel; ++i)
467  {
468  wk1[i] = Array<OneD, NekDouble> (physTot, 0.0);
469  wk2[i] = Array<OneD, NekDouble> (physTot, 0.0);
470  gradP[i] = Array<OneD, NekDouble> (physTot, 0.0);
471  }
472 
473  // Jacobian
474  Array<OneD, NekDouble> Jac(physTot, 0.0);
475  m_mapping->GetJacobian(Jac);
476 
477  // Factors for Laplacian system
479  factors[StdRegions::eFactorLambda] = 0.0;
480 
481  m_pressure->BwdTrans(m_pressure->GetCoeffs(),
482  m_pressure->UpdatePhys());
483  forcing_L2 = m_pressure->L2(Forcing, wk1[0]);
484  while (!converged)
485  {
486  // Update iteration counter and set previous iteration field
487  // (use previous timestep solution for first iteration)
488  s++;
489  ASSERTL0(s < maxIter,
490  "VCSMapping exceeded maximum number of iterations.");
491 
492  Vmath::Vcopy(physTot, m_pressure->GetPhys(), 1,
493  previous_iter, 1);
494 
495  // Correct pressure bc to account for iteration
496  m_extrapolation->CorrectPressureBCs(previous_iter);
497 
498  //
499  // Calculate forcing term for this iteration
500  //
501  for(int i = 0; i < nvel; ++i)
502  {
504  previous_iter, gradP[i]);
505  if(m_pressure->GetWaveSpace())
506  {
507  m_pressure->HomogeneousBwdTrans(gradP[i], wk1[i]);
508  }
509  else
510  {
511  Vmath::Vcopy(physTot, gradP[i], 1, wk1[i], 1);
512  }
513  }
514  m_mapping->RaiseIndex(wk1, wk2); // G(p)
515 
516  m_mapping->Divergence(wk2, F_corrected); // div(G(p))
517  if (!m_mapping->HasConstantJacobian())
518  {
519  Vmath::Vmul(physTot, F_corrected, 1,
520  Jac, 1,
521  F_corrected, 1);
522  }
523  // alpha*J*div(G(p))
524  Vmath::Smul(physTot, m_pressureRelaxation, F_corrected, 1,
525  F_corrected, 1);
526  if(m_pressure->GetWaveSpace())
527  {
528  m_pressure->HomogeneousFwdTrans(F_corrected, F_corrected);
529  }
530  // alpha*J*div(G(p)) - p_ii
531  for (int i = 0; i < m_nConvectiveFields; ++i)
532  {
534  gradP[i], wk1[0]);
535  Vmath::Vsub(physTot, F_corrected, 1, wk1[0], 1,
536  F_corrected, 1);
537  }
538  // p_i,i - J*div(G(p))
539  Vmath::Neg(physTot, F_corrected, 1);
540  // alpha*F - alpha*J*div(G(p)) + p_i,i
541  Vmath::Smul(physTot, m_pressureRelaxation, Forcing, 1,
542  wk1[0], 1);
543  Vmath::Vadd(physTot, wk1[0], 1, F_corrected, 1, F_corrected, 1);
544 
545  //
546  // Solve system
547  //
548  m_pressure->HelmSolve(F_corrected, m_pressure->UpdateCoeffs(),
549  NullFlagList, factors);
550  m_pressure->BwdTrans(m_pressure->GetCoeffs(),
551  m_pressure->UpdatePhys());
552 
553  //
554  // Test convergence
555  //
556  error = m_pressure->L2(m_pressure->GetPhys(), previous_iter);
557  if ( forcing_L2 != 0)
558  {
559  if ( (error/forcing_L2 < m_pressureTolerance))
560  {
561  converged = true;
562  }
563  }
564  else
565  {
566  if ( error < m_pressureTolerance)
567  {
568  converged = true;
569  }
570  }
571  }
572  if (m_verbose && m_session->GetComm()->GetRank()==0)
573  {
574  std::cout << " Pressure system (mapping) converged in " << s <<
575  " iterations with error = " << error << std::endl;
576  }
577  }
578  }
579 
580  /**
581  * Solve velocity system
582  */
584  const Array<OneD, const Array<OneD, NekDouble> > &Forcing,
585  Array<OneD, Array<OneD, NekDouble> > &outarray,
586  const NekDouble aii_Dt)
587  {
588  if(!m_implicitViscous)
589  {
591  }
592  else
593  {
594  int physTot = m_fields[0]->GetTotPoints();
595  int nvel = m_nConvectiveFields;
596  bool converged = false; // flag to mark if system converged
597  int s = 0; // iteration counter
598  NekDouble error, max_error; // L2 error at current iteration
599 
600  int maxIter;
601  m_session->LoadParameter("MappingMaxIter",maxIter,5000);
602 
603  //L2 norm of F
605 
606  // rhs of the equation at current iteration
607  Array<OneD, Array<OneD, NekDouble> > F_corrected(nvel);
608  // Solution at previous iteration
609  Array<OneD, Array<OneD, NekDouble> > previous_iter(nvel);
610  // Working space
612  for(int i = 0; i < nvel; ++i)
613  {
614  F_corrected[i] = Array<OneD, NekDouble> (physTot, 0.0);
615  previous_iter[i] = Array<OneD, NekDouble> (physTot, 0.0);
616  wk[i] = Array<OneD, NekDouble> (physTot, 0.0);
617  }
618 
619  // Factors for Helmholtz system
621  factors[StdRegions::eFactorLambda] =
622  1.0*m_viscousRelaxation/aii_Dt/m_kinvis;
623  if(m_useSpecVanVisc)
624  {
628  }
629 
630  // Calculate L2-norm of F and set initial solution for iteration
631  for(int i = 0; i < nvel; ++i)
632  {
633  forcing_L2[i] = m_fields[0]->L2(Forcing[i],wk[0]);
634  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
635  previous_iter[i]);
636  }
637 
638  while (!converged)
639  {
640  converged = true;
641  // Iteration counter
642  s++;
643  ASSERTL0(s < maxIter,
644  "VCSMapping exceeded maximum number of iterations.");
645 
646  max_error = 0.0;
647 
648  //
649  // Calculate forcing term for next iteration
650  //
651 
652  // Calculate L(U)- in this parts all components might be coupled
653  if(m_fields[0]->GetWaveSpace())
654  {
655  for (int i = 0; i < nvel; ++i)
656  {
657  m_fields[0]->HomogeneousBwdTrans(previous_iter[i],
658  wk[i]);
659  }
660  }
661  else
662  {
663  for (int i = 0; i < nvel; ++i)
664  {
665  Vmath::Vcopy(physTot, previous_iter[i], 1, wk[i], 1);
666  }
667  }
668 
669  // (L(U^i) - 1/alpha*U^i_jj)
670  m_mapping->VelocityLaplacian(wk, F_corrected,
671  1.0/m_viscousRelaxation);
672 
673  if(m_fields[0]->GetWaveSpace())
674  {
675  for (int i = 0; i < nvel; ++i)
676  {
677  m_fields[0]->HomogeneousFwdTrans(F_corrected[i],
678  F_corrected[i]);
679  }
680  }
681  else
682  {
683  for (int i = 0; i < nvel; ++i)
684  {
685  Vmath::Vcopy(physTot, F_corrected[i], 1,
686  F_corrected[i], 1);
687  }
688  }
689 
690  // Loop velocity components
691  for (int i = 0; i < nvel; ++i)
692  {
693  // (-alpha*L(U^i) + U^i_jj)
694  Vmath::Smul(physTot, -1.0*m_viscousRelaxation,
695  F_corrected[i], 1,
696  F_corrected[i], 1);
697  // F_corrected = alpha*F + (-alpha*L(U^i) + U^i_jj)
698  Vmath::Smul(physTot, m_viscousRelaxation, Forcing[i], 1,
699  wk[0], 1);
700  Vmath::Vadd(physTot, wk[0], 1, F_corrected[i], 1,
701  F_corrected[i], 1);
702 
703  //
704  // Solve System
705  //
706  m_fields[i]->HelmSolve(F_corrected[i],
707  m_fields[i]->UpdateCoeffs(),
708  NullFlagList, factors);
709  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),outarray[i]);
710 
711  //
712  // Test convergence
713  //
714  error = m_fields[i]->L2(outarray[i], previous_iter[i]);
715 
716  if ( forcing_L2[i] != 0)
717  {
718  if ( (error/forcing_L2[i] >= m_viscousTolerance))
719  {
720  converged = false;
721  }
722  }
723  else
724  {
725  if ( error >= m_viscousTolerance)
726  {
727  converged = false;
728  }
729  }
730  if (error > max_error)
731  {
732  max_error = error;
733  }
734 
735  // Copy field to previous_iter
736  Vmath::Vcopy(physTot, outarray[i], 1, previous_iter[i], 1);
737  }
738  }
739  if (m_verbose && m_session->GetComm()->GetRank()==0)
740  {
741  std::cout << " Velocity system (mapping) converged in " << s <<
742  " iterations with error = " << max_error << std::endl;
743  }
744  }
745  }
746 
747  /**
748  * Explicit terms of the mapping
749  */
751  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
752  Array<OneD, Array<OneD, NekDouble> > &outarray)
753  {
754  int physTot = m_fields[0]->GetTotPoints();
759  for (int i = 0; i < m_nConvectiveFields; ++i)
760  {
761  velPhys[i] = Array<OneD, NekDouble> (physTot, 0.0);
762  Forcing[i] = Array<OneD, NekDouble> (physTot, 0.0);
763  tmp[i] = Array<OneD, NekDouble> (physTot, 0.0);
764  }
765 
766  // Get fields and store velocity in wavespace and physical space
767  if(m_fields[0]->GetWaveSpace())
768  {
769  for (int i = 0; i < m_nConvectiveFields; ++i)
770  {
771  vel[i] = inarray[i];
772  m_fields[0]->HomogeneousBwdTrans(vel[i],velPhys[i]);
773  }
774  }
775  else
776  {
777  for (int i = 0; i < m_nConvectiveFields; ++i)
778  {
779  vel[i] = inarray[i];
780  Vmath::Vcopy(physTot, inarray[i], 1, velPhys[i], 1);
781  }
782  }
783 
784  //Advection contribution
785  MappingAdvectionCorrection(velPhys, Forcing);
786 
787  // Time-derivative contribution
788  if ( m_mapping->IsTimeDependent() )
789  {
790  MappingAccelerationCorrection(vel, velPhys, tmp);
791  for (int i = 0; i < m_nConvectiveFields; ++i)
792  {
793  Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
794  }
795  }
796 
797  // Pressure contribution
798  if (!m_implicitPressure)
799  {
801  for (int i = 0; i < m_nConvectiveFields; ++i)
802  {
803  Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
804  }
805  }
806  // Viscous contribution
807  if ( (!m_implicitViscous) && (!m_neglectViscous))
808  {
809  MappingViscousCorrection(velPhys, tmp);
810  for (int i = 0; i < m_nConvectiveFields; ++i)
811  {
812  Vmath::Smul(physTot, m_kinvis, tmp[i], 1, tmp[i], 1);
813  Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
814  }
815  }
816 
817  // If necessary, transform to wavespace
818  if(m_fields[0]->GetWaveSpace())
819  {
820  for (int i = 0; i < m_nConvectiveFields; ++i)
821  {
822  m_fields[0]->HomogeneousFwdTrans(Forcing[i],Forcing[i]);
823  }
824  }
825 
826  // Add to outarray
827  for (int i = 0; i < m_nConvectiveFields; ++i)
828  {
829  Vmath::Vadd(physTot, outarray[i], 1, Forcing[i], 1, outarray[i], 1);
830  }
831  }
832 
834  const Array<OneD, const Array<OneD, NekDouble> > &velPhys,
835  Array<OneD, Array<OneD, NekDouble> > &outarray)
836  {
837  int physTot = m_fields[0]->GetTotPoints();
838  int nvel = m_nConvectiveFields;
839 
840  Array<OneD, Array<OneD, NekDouble> > wk(nvel*nvel);
841 
842  // Apply Christoffel symbols to obtain {i,kj}vel(k)
843  m_mapping->ApplyChristoffelContravar(velPhys, wk);
844 
845  // Calculate correction -U^j*{i,kj}vel(k)
846  for (int i = 0; i< nvel; i++)
847  {
848  Vmath::Zero(physTot,outarray[i],1);
849  for (int j = 0; j< nvel; j++)
850  {
851  Vmath::Vvtvp(physTot,wk[i*nvel+j],1,velPhys[j],1,
852  outarray[i],1,outarray[i],1);
853  }
854  Vmath::Neg(physTot, outarray[i], 1);
855  }
856  }
857 
859  const Array<OneD, const Array<OneD, NekDouble> > &vel,
860  const Array<OneD, const Array<OneD, NekDouble> > &velPhys,
861  Array<OneD, Array<OneD, NekDouble> > &outarray)
862  {
863  int physTot = m_fields[0]->GetTotPoints();
864  int nvel = m_nConvectiveFields;
865 
866  Array<OneD, Array<OneD, NekDouble> > wk(nvel*nvel);
868  Array<OneD, Array<OneD, NekDouble> > coordVel(nvel);
869  for (int i = 0; i< nvel; i++)
870  {
871  tmp[i] = Array<OneD, NekDouble> (physTot, 0.0);
872  coordVel[i] = Array<OneD, NekDouble> (physTot, 0.0);
873  }
874  // Get coordinates velocity in transformed system
875  m_mapping->GetCoordVelocity(tmp);
876  m_mapping->ContravarFromCartesian(tmp, coordVel);
877 
878  // Calculate first term: U^j u^i,j = U^j (du^i/dx^j + {i,kj}u^k)
879  m_mapping->ApplyChristoffelContravar(velPhys, wk);
880  for (int i=0; i< nvel; i++)
881  {
882  Vmath::Zero(physTot,outarray[i],1);
883 
884  m_fields[0]->PhysDeriv(velPhys[i], tmp[0], tmp[1]);
885  for (int j=0; j< nvel; j++)
886  {
887  if (j == 2)
888  {
889  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j],
890  vel[i], tmp[2]);
891  if (m_fields[0]->GetWaveSpace())
892  {
893  m_fields[0]->HomogeneousBwdTrans(tmp[2],tmp[2]);
894  }
895  }
896 
897  Vmath::Vadd(physTot,wk[i*nvel+j],1,tmp[j],1,
898  wk[i*nvel+j], 1);
899 
900  Vmath::Vvtvp(physTot, coordVel[j], 1, wk[i*nvel+j], 1,
901  outarray[i], 1, outarray[i], 1);
902  }
903  }
904 
905  // Set wavespace to false and store current value
906  bool wavespace = m_fields[0]->GetWaveSpace();
907  m_fields[0]->SetWaveSpace(false);
908 
909  // Add -u^j U^i,j
910  m_mapping->ApplyChristoffelContravar(coordVel, wk);
911  for (int i=0; i< nvel; i++)
912  {
913  if(nvel == 2)
914  {
915  m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1]);
916  }
917  else
918  {
919  m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1], tmp[2]);
920  }
921 
922  for (int j=0; j< nvel; j++)
923  {
924  Vmath::Vadd(physTot,wk[i*nvel+j],1,tmp[j],1,
925  wk[i*nvel+j], 1);
926  Vmath::Neg(physTot, wk[i*nvel+j], 1);
927 
928  Vmath::Vvtvp(physTot, velPhys[j], 1, wk[i*nvel+j], 1,
929  outarray[i], 1, outarray[i], 1);
930  }
931  }
932 
933  // Restore value of wavespace
934  m_fields[0]->SetWaveSpace(wavespace);
935  }
936 
938  Array<OneD, Array<OneD, NekDouble> > &outarray)
939  {
940  int physTot = m_fields[0]->GetTotPoints();
941  int nvel = m_nConvectiveFields;
942 
943  // Calculate g^(ij)p_(,j)
944  m_mapping->RaiseIndex(m_gradP, outarray);
945 
946  // Calculate correction = (nabla p)/J - g^(ij)p_,j
947  // (Jac is not required if it is constant)
948  if ( !m_mapping->HasConstantJacobian())
949  {
950  Array<OneD, NekDouble> Jac(physTot, 0.0);
951  m_mapping->GetJacobian(Jac);
952  for(int i = 0; i < nvel; ++i)
953  {
954  Vmath::Vdiv(physTot, m_gradP[i], 1, Jac, 1, m_gradP[i], 1);
955  }
956  }
957  for(int i = 0; i < nvel; ++i)
958  {
959  Vmath::Vsub(physTot, m_gradP[i], 1,outarray[i], 1,
960  outarray[i],1);
961  }
962  }
963 
965  const Array<OneD, const Array<OneD, NekDouble> > &velPhys,
966  Array<OneD, Array<OneD, NekDouble> > &outarray)
967  {
968  // L(U) - 1.0*d^2(u^i)/dx^jdx^j
969  m_mapping->VelocityLaplacian(velPhys, outarray, 1.0);
970  }
971 
972 } //end of namespace
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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:750
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:858
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:583
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:833
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:442
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:362
STL namespace.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
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
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:964
IMEX 2nd order scheme using Backward Different Formula & Extrapolation.
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
Definition: VCSMapping.cpp:439
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:236
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
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:396
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()
Array< OneD, Array< OneD, NekDouble > > m_F
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
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:937
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:373
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:1061
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
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: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
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:190
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