Nektar++
MMFAdvection.cpp
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////
2 //
3 // File MMFAdvection.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: MMF solve routines
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <iomanip>
36 #include <iostream>
37 
38 #include <boost/algorithm/string.hpp>
39 #include <boost/core/ignore_unused.hpp>
40 
44 
45 namespace Nektar
46 {
47 std::string MMFAdvection::className =
49  "MMFAdvection", MMFAdvection::create, "MMFAdvection equation.");
50 
53  : UnsteadySystem(pSession, pGraph), MMFSystem(pSession, pGraph),
54  AdvectionSystem(pSession, pGraph)
55 {
56  m_planeNumber = 0;
57 }
58 
59 /**
60  * @brief Initialisation object for the unsteady linear advection equation.
61  */
62 void MMFAdvection::v_InitObject(bool DeclareFields)
63 {
64  // Call to the initialisation object
65  UnsteadySystem::v_InitObject(DeclareFields);
66 
67  int nq = m_fields[0]->GetNpoints();
68  int shapedim = m_fields[0]->GetShapeDimension();
69  Array<OneD, Array<OneD, NekDouble>> Anisotropy(shapedim);
70  for (int j = 0; j < shapedim; ++j)
71  {
72  Anisotropy[j] = Array<OneD, NekDouble>(nq, 1.0);
73  }
74 
75  MMFSystem::MMFInitObject(Anisotropy);
76 
77  // Define TestType
78  ASSERTL0(m_session->DefinesSolverInfo("TESTTYPE"),
79  "No TESTTYPE defined in session.");
80  std::string TestTypeStr = m_session->GetSolverInfo("TESTTYPE");
81  for (int i = 0; i < (int)SIZE_TestType; ++i)
82  {
83  if (boost::iequals(TestTypeMap[i], TestTypeStr))
84  {
85  m_TestType = (TestType)i;
86  break;
87  }
88  }
89 
90  m_session->LoadParameter("Angular Frequency", m_waveFreq, m_pi);
91  m_session->LoadParameter("Rotational Angle", m_RotAngle, 0.0);
92  m_session->LoadParameter("Velocity Projection", m_VelProjection, 0);
93 
94  // Read the advection velocities from session file
95  m_session->LoadParameter("advx", m_advx, 1.0);
96  m_session->LoadParameter("advy", m_advy, 1.0);
97  m_session->LoadParameter("advz", m_advz, 1.0);
98 
99  std::vector<std::string> vel;
100  vel.push_back("Vx");
101  vel.push_back("Vy");
102  vel.push_back("Vz");
103 
104  // Resize the advection velocities vector to dimension of the problem
105  vel.resize(m_spacedim);
106 
107  // Store in the global variable m_velocity the advection velocities
108 
110  for (int k = 0; k < m_spacedim; ++k)
111  {
113  }
114 
115  switch (m_surfaceType)
116  {
121  {
122  // true = project velocity onto the tangent plane
124  }
125  break;
126 
127  case SolverUtils::ePlane:
128  case SolverUtils::eCube:
129  {
130  GetFunction("AdvectionVelocity")->Evaluate(vel, m_velocity);
131  }
132  break;
133 
134  default:
135  break;
136  }
137 
138  std::cout << "|Velocity vector| = ( " << RootMeanSquare(m_velocity[0])
139  << " , " << RootMeanSquare(m_velocity[1]) << " , "
140  << RootMeanSquare(m_velocity[2]) << " ) " << std::endl;
141 
142  // Define the normal velocity fields
143  if (m_fields[0]->GetTrace())
144  {
146  }
147 
148  std::string advName;
149  std::string riemName;
150  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
151  m_advObject =
153  m_advObject->SetFluxVector(&MMFAdvection::GetFluxVector, this);
154  m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
156  riemName, m_session);
157  m_riemannSolver->SetScalar("Vn", &MMFAdvection::GetNormalVelocity, this);
158 
159  m_advObject->SetRiemannSolver(m_riemannSolver);
160  m_advObject->InitObject(m_session, m_fields);
161 
162  // Compute m_traceVn = n \cdot v
164 
165  // Compute m_vellc = nabal a^j \cdot m_vel
167  std::cout << "m_vellc is generated with mag = " << AvgInt(m_vellc)
168  << std::endl;
169 
170  // Compute vel \cdot MF
172 
173  // Modify e^i as v^i e^i
174  for (int j = 0; j < m_shapedim; j++)
175  {
176  for (int k = 0; k < m_spacedim; k++)
177  {
178  Vmath::Vmul(nq, &m_veldotMF[j][0], 1, &m_movingframes[j][k * nq], 1,
179  &m_movingframes[j][k * nq], 1);
180  }
181  }
182 
183  // Reflect it into m_ncdotMFFwd and Bwd
184  ComputencdotMF();
185 
186  // If explicit it computes RHS and PROJECTION for the time integration
188  {
191  }
192  // Otherwise it gives an error (no implicit integration)
193  else
194  {
195  ASSERTL0(false, "Implicit unsteady Advection not set up.");
196  }
197 }
198 
199 /**
200  * @brief Unsteady linear advection equation destructor.
201  */
203 {
204 }
205 
207 {
208  ASSERTL0(m_intScheme != 0, "No time integration scheme.");
209 
210  int i, nchk = 1;
211  int nvariables = 0;
212  int nfields = m_fields.size();
213  int nq = m_fields[0]->GetNpoints();
214 
215  if (m_intVariables.empty())
216  {
217  for (i = 0; i < nfields; ++i)
218  {
219  m_intVariables.push_back(i);
220  }
221  nvariables = nfields;
222  }
223  else
224  {
225  nvariables = m_intVariables.size();
226  }
227 
228  // Set up wrapper to fields data storage.
229  Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
230  Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
231 
232  // Order storage to list time-integrated fields first.
233  for (i = 0; i < nvariables; ++i)
234  {
235  fields[i] = m_fields[m_intVariables[i]]->GetPhys();
236  m_fields[m_intVariables[i]]->SetPhysState(false);
237  }
238 
239  // Initialise time integration scheme
240  m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
241 
242  // Check uniqueness of checkpoint output
243  ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
244  (m_checktime > 0.0 && m_checksteps == 0) ||
245  (m_checktime == 0.0 && m_checksteps > 0),
246  "Only one of IO_CheckTime and IO_CheckSteps "
247  "should be set!");
248 
249  LibUtilities::Timer timer;
250  bool doCheckTime = false;
251  int step = 0;
252  NekDouble intTime = 0.0;
253  NekDouble cpuTime = 0.0;
254  NekDouble elapsed = 0.0;
255 
256  int Ntot, indx;
257  // Perform integration in time.
258  Ntot = m_steps / m_checksteps + 1;
259 
260  Array<OneD, NekDouble> dMass(Ntot);
261 
262  Array<OneD, NekDouble> zeta(nq);
263  Array<OneD, Array<OneD, NekDouble>> fieldsprimitive(nvariables);
264  for (int i = 0; i < nvariables; ++i)
265  {
266  fieldsprimitive[i] = Array<OneD, NekDouble>(nq);
267  }
268 
269  while (step < m_steps || m_time < m_fintime - NekConstants::kNekZeroTol)
270  {
271  timer.Start();
272  fields = m_intScheme->TimeIntegrate(step, m_timestep, m_ode);
273  timer.Stop();
274 
275  m_time += m_timestep;
276  elapsed = timer.TimePerTest(1);
277  intTime += elapsed;
278  cpuTime += elapsed;
279 
280  // Write out status information
281  if (m_session->GetComm()->GetRank() == 0 && !((step + 1) % m_infosteps))
282  {
283  std::cout << "Steps: " << std::setw(8) << std::left << step + 1
284  << " "
285  << "Time: " << std::setw(12) << std::left << m_time;
286 
287  std::stringstream ss;
288  ss << cpuTime << "s";
289  std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
290  << std::endl;
291 
292  // Masss = h^*
293  indx = (step + 1) / m_checksteps;
294  dMass[indx] =
295  (m_fields[0]->Integral(fields[0]) - m_Mass0) / m_Mass0;
296 
297  std::cout << "dMass = " << std::setw(8) << std::left << dMass[indx]
298  << std::endl;
299 
300  cpuTime = 0.0;
301  }
302 
303  // Transform data into coefficient space
304  for (i = 0; i < nvariables; ++i)
305  {
306  m_fields[m_intVariables[i]]->SetPhys(fields[i]);
307  m_fields[m_intVariables[i]]->FwdTransLocalElmt(
308  fields[i], m_fields[m_intVariables[i]]->UpdateCoeffs());
309  m_fields[m_intVariables[i]]->SetPhysState(false);
310  }
311 
312  // Write out checkpoint files
313  if ((m_checksteps && step && !((step + 1) % m_checksteps)) ||
314  doCheckTime)
315  {
316  Checkpoint_Output(nchk++);
317  doCheckTime = false;
318  }
319 
320  // Step advance
321  ++step;
322  }
323 
324  std::cout << "dMass = ";
325  for (i = 0; i < Ntot; ++i)
326  {
327  std::cout << dMass[i] << " , ";
328  }
329  std::cout << std::endl << std::endl;
330 
331  // Print out summary statistics
332  if (m_session->GetComm()->GetRank() == 0)
333  {
334  if (m_cflSafetyFactor > 0.0)
335  {
336  std::cout << "CFL safety factor : " << m_cflSafetyFactor
337  << std::endl
338  << "CFL time-step : " << m_timestep << std::endl;
339  }
340 
341  if (m_session->GetSolverInfo("Driver") != "SteadyState")
342  {
343  std::cout << "Time-integration : " << intTime << "s" << std::endl;
344  }
345  }
346 
347  for (i = 0; i < nvariables; ++i)
348  {
349  m_fields[m_intVariables[i]]->SetPhys(fields[i]);
350  m_fields[m_intVariables[i]]->SetPhysState(true);
351  }
352 
353  for (i = 0; i < nvariables; ++i)
354  {
355  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
356  m_fields[i]->UpdateCoeffs());
357  }
358 }
359 
360 /**
361  * @brief Get the normal velocity for the linear advection equation.
362  */
364 {
365  // Number of trace (interface) points
366  int i;
367  int nTracePts = GetTraceNpoints();
368 
369  // Auxiliary variable to compute the normal velocity
370  Array<OneD, NekDouble> tmp(nTracePts);
371 
372  // Reset the normal velocity
373  Vmath::Zero(nTracePts, m_traceVn, 1);
374 
375  for (i = 0; i < m_velocity.size(); ++i)
376  {
377  m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
378 
379  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, tmp, 1, m_traceVn, 1,
380  m_traceVn, 1);
381  }
382 
383  return m_traceVn;
384 }
385 
386 /**
387  * @brief Compute the right-hand side for the linear advection equation.
388  *
389  * @param inarray Given fields.
390  * @param outarray Calculated solution.
391  * @param time Time.
392  */
394  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
395  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
396 {
397  boost::ignore_unused(time);
398 
399  int i;
400  int nvariables = inarray.size();
401  int npoints = GetNpoints();
402 
403  switch (m_projectionType)
404  {
406  {
407  int ncoeffs = inarray[0].size();
408 
409  if (m_spacedim == 3)
410  {
411  Array<OneD, Array<OneD, NekDouble>> WeakAdv(nvariables);
412 
413  WeakAdv[0] = Array<OneD, NekDouble>(ncoeffs * nvariables);
414  for (i = 1; i < nvariables; ++i)
415  {
416  WeakAdv[i] = WeakAdv[i - 1] + ncoeffs;
417  }
418 
419  // Compute \nabla \cdot \vel u according to MMF scheme
420  WeakDGDirectionalAdvection(inarray, WeakAdv);
421 
422  for (i = 0; i < nvariables; ++i)
423  {
424  m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
425  m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
426 
427  // Add m_vellc * inarray[i] = \nabla v^m \cdot e^m to
428  // outarray[i]
429  // Vmath::Vvtvp(npoints, &m_vellc[0], 1, &inarray[i][0], 1,
430  // &outarray[i][0], 1, &outarray[i][0], 1);
431  Vmath::Neg(npoints, outarray[i], 1);
432  }
433  }
434  else
435  {
436  m_advObject->Advect(2, m_fields, m_velocity, inarray, outarray,
437  0.0);
438 
439  for (i = 0; i < nvariables; ++i)
440  {
441  Vmath::Neg(npoints, outarray[i], 1);
442  }
443  }
444  }
445  break;
446 
447  default:
448  {
449  ASSERTL0(false, "Unknown projection scheme");
450  }
451  break;
452  }
453 }
454 /**
455  * @brief Compute the projection for the linear advection equation.
456  *
457  * @param inarray Given fields.
458  * @param outarray Calculated solution.
459  * @param time Time.
460  */
462  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
463  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
464 {
465  // Counter variable
466  int i;
467 
468  // Number of fields (variables of the problem)
469  int nVariables = inarray.size();
470 
471  // Set the boundary conditions
472  SetBoundaryConditions(time);
473 
474  // Switch on the projection type (Discontinuous or Continuous)
475  switch (m_projectionType)
476  {
477  // Discontinuous projection
479  {
480  // Number of quadrature points
481  int nQuadraturePts = GetNpoints();
482 
483  // Just copy over array
484  for (i = 0; i < nVariables; ++i)
485  {
486  Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
487  }
488  break;
489  }
490 
491  // Continuous projection
494  {
495  Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs(), 0.0);
496  for (i = 0; i < nVariables; ++i)
497  {
498  m_fields[i]->FwdTrans(inarray[i], coeffs);
499  m_fields[i]->BwdTrans(coeffs, outarray[i]);
500  }
501  break;
502  }
503 
504  default:
505  ASSERTL0(false, "Unknown projection scheme");
506  break;
507  }
508 }
509 
510 /**
511  * @brief Return the flux vector for the linear advection equation.
512  *
513  * @param i Component of the flux vector to calculate.
514  * @param physfield Fields.
515  * @param flux Resulting flux.
516  */
518  const Array<OneD, Array<OneD, NekDouble>> &physfield,
520 {
521  ASSERTL1(flux[0].size() == m_velocity.size(),
522  "Dimension of flux array and velocity array do not match");
523 
524  int i, j;
525  int nq = physfield[0].size();
526 
527  for (i = 0; i < flux.size(); ++i)
528  {
529  for (j = 0; j < flux[0].size(); ++j)
530  {
531  Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1, flux[i][j], 1);
532  }
533  }
534 }
535 
537  const Array<OneD, Array<OneD, NekDouble>> &InField,
538  Array<OneD, Array<OneD, NekDouble>> &OutField)
539 {
540  int i, j;
541  int ncoeffs = GetNcoeffs();
542  int nTracePointsTot = GetTraceNpoints();
543  int nvariables = m_fields.size();
544 
545  Array<OneD, Array<OneD, NekDouble>> physfield(nvariables);
546 
547  // Get the variables in physical space
548  // already in physical space
549  for (i = 0; i < nvariables; ++i)
550  {
551  physfield[i] = InField[i];
552  }
553 
555  for (i = 0; i < nvariables; ++i)
556  {
557  for (j = 0; j < m_shapedim; ++j)
558  {
559  WeakDeriv[j] = Array<OneD, NekDouble>(ncoeffs, 0.0);
560 
561  // Directional derivation with respect to the j'th moving frame
562  // tmp[j] = \nabla \physfield[i] \cdot \mathbf{e}^j
563  // Implemented at TriExp::v_IProductWRTDirectionalDerivBase_SumFa
564  m_fields[0]->IProductWRTDirectionalDerivBase(
565  m_movingframes[j], physfield[i], WeakDeriv[j]);
566  }
567 
568  // Get the numerical flux and add to the modal coeffs
569  // if the NumericalFluxs function already includes the
570  // normal in the output
571 
572  Array<OneD, NekDouble> Fwd(nTracePointsTot);
573  Array<OneD, NekDouble> Bwd(nTracePointsTot);
574 
575  Array<OneD, NekDouble> flux(nTracePointsTot, 0.0);
576  Array<OneD, NekDouble> fluxFwd(nTracePointsTot);
577  Array<OneD, NekDouble> fluxBwd(nTracePointsTot);
578 
579  // Evaluate numerical flux in physical space which may in
580  // general couple all component of vectors
581  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
582 
583  // evaulate upwinded m_fields[i]
584  m_fields[i]->GetTrace()->Upwind(m_traceVn, Fwd, Bwd, flux);
585 
586  OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
587  for (j = 0; j < m_shapedim; ++j)
588  {
589  // calculate numflux = (n \cdot MF)*flux
590  Vmath::Vmul(nTracePointsTot, &flux[0], 1, &m_ncdotMFFwd[j][0], 1,
591  &fluxFwd[0], 1);
592  Vmath::Vmul(nTracePointsTot, &flux[0], 1, &m_ncdotMFBwd[j][0], 1,
593  &fluxBwd[0], 1);
594  Vmath::Neg(ncoeffs, WeakDeriv[j], 1);
595 
596  // FwdBwdtegral because generallize (N \cdot MF)_{FWD} \neq -(N
597  // \cdot MF)_{BWD}
598  m_fields[i]->AddFwdBwdTraceIntegral(fluxFwd, fluxBwd, WeakDeriv[j]);
599  m_fields[i]->SetPhysState(false);
600 
601  Vmath::Vadd(ncoeffs, &WeakDeriv[j][0], 1, &OutField[i][0], 1,
602  &OutField[i][0], 1);
603  }
604  }
605 }
606 
608  Array<OneD, Array<OneD, NekDouble>> &velocity)
609 {
610  int nq = m_fields[0]->GetNpoints();
611 
612  NekDouble vel_phi, vel_theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
613  NekDouble x0j, x1j, x2j;
614 
615  Array<OneD, NekDouble> x0(nq);
616  Array<OneD, NekDouble> x1(nq);
617  Array<OneD, NekDouble> x2(nq);
618 
619  m_fields[0]->GetCoords(x0, x1, x2);
620 
621  // theta = a*sin(z/r), phi = a*tan(y/x);
622  for (int j = 0; j < nq; j++)
623  {
624  x0j = x0[j];
625  x1j = x1[j];
626  x2j = x2[j];
627 
628  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
629  cos_theta);
630 
631  vel_phi = m_waveFreq * (cos_theta * cos(m_RotAngle) +
632  sin_theta * cos_varphi * sin(m_RotAngle));
633  vel_theta = -1.0 * m_waveFreq * sin_theta * sin(m_RotAngle);
634 
635  velocity[0][j] =
636  -1.0 * vel_phi * sin_varphi - vel_theta * sin_theta * cos_varphi;
637  velocity[1][j] =
638  vel_phi * cos_varphi - vel_theta * sin_theta * sin_varphi;
639  velocity[2][j] = vel_theta * cos_theta;
640  }
641 
642  // Project the veloicty on the tangent plane
643 
644  if (m_VelProjection)
645  {
646  // Check MovingFrames \cdot SurfaceNormals = 0
648 
652 
653  for (int k = 0; k < m_spacedim; ++k)
654  {
655  newvelocity[k] = Array<OneD, NekDouble>(nq);
656  MF1[k] = Array<OneD, NekDouble>(nq);
657  MF2[k] = Array<OneD, NekDouble>(nq);
658  SN[k] = Array<OneD, NekDouble>(nq);
659 
660  Vmath::Vcopy(nq, &m_movingframes[0][k * nq], 1, &MF1[k][0], 1);
661  Vmath::Vcopy(nq, &m_movingframes[1][k * nq], 1, &MF2[k][0], 1);
662  }
663 
664  VectorCrossProd(MF1, MF2, SN);
665  GramSchumitz(SN, m_velocity, newvelocity, true);
666 
667  Array<OneD, NekDouble> tmp(nq, 0.0);
669 
670  for (int k = 0; k < m_spacedim; ++k)
671  {
672  Vmath::Vsub(nq, &m_velocity[0][0], 1, &newvelocity[0][0], 1,
673  &tmp[0], 1);
674  Vmath::Vmul(nq, &tmp[0], 1, &tmp[0], 1, &tmp[0], 1);
675  Vmath::Vadd(nq, tmp, 1, Projection, 1, Projection, 1);
676  }
677 
678  std::cout
679  << "Velocity vector is projected onto the tangent plane: Diff = "
680  << RootMeanSquare(Projection) << std::endl;
681 
682  for (int k = 0; k < m_spacedim; ++k)
683  {
684  Vmath::Vcopy(nq, &newvelocity[k][0], 1, &m_velocity[k][0], 1);
685  }
686  }
687 }
688 
689 /*
690 void MMFAdvection::EvaluateAdvectionVelocity(Array<OneD, Array<OneD, NekDouble>
691 > &velocity)
692  {
693  int nq = m_fields[0]->GetNpoints();
694 
695  NekDouble vel_phi, sin_phi, cos_phi;
696  NekDouble vel_theta, sin_theta, cos_theta;
697  NekDouble radius, xyrad, Tol = 0.00001;
698 
699  Array<OneD,NekDouble> x0(nq);
700  Array<OneD,NekDouble> x1(nq);
701  Array<OneD,NekDouble> x2(nq);
702 
703  m_fields[0]->GetCoords(x0,x1,x2);
704 
705  // theta = a*sin(z/r), phi = a*tan(y/x);
706  for (int j = 0; j < nq; j++)
707  {
708  switch(m_surfaceType)
709  {
710  case SolverUtils::eSphere:
711  case SolverUtils::eTRSphere:
712  {
713  radius = sqrt( x0[j]*x0[j] + x1[j]*x1[j] + x2[j]*x2[j] );
714  }
715  break;
716 
717  case SolverUtils::eIrregular:
718  {
719  radius = sqrt( 2.0*x0[j]*x0[j] + x1[j]*x1[j]*x1[j]*x1[j] +
720 x1[j]*x1[j]
721  + x2[j]*x2[j]*x2[j]*x2[j] + x2[j]*x2[j] );
722  }
723  break;
724 
725  // 2 x^2 + 2(y^4 - y ) + z^4 + z^2 = 2.0
726  case SolverUtils::eNonconvex:
727  {
728  radius = sqrt( 2.0*x0[j]*x0[j] + 2.0*( x1[j]*x1[j]*x1[j]*x1[j] -
729 x1[j]*x1[j] )
730  + x2[j]*x2[j]*x2[j]*x2[j] + x2[j]*x2[j] );
731  }
732  break;
733 
734  default:
735  break;
736  }
737 
738  // At North and South poles, (ax,ay,ax) = (0, Omega_f*sin(alpha),0)
739  if( fabs(fabs(x2[j]) - radius) < Tol )
740  {
741  sin_theta = x2[j]/radius;
742 
743  velocity[0][j] = 0.0;
744  velocity[1][j] = 0.0;
745  velocity[2][j] = 0.0;
746  }
747  else
748  {
749  // Compute the arc length of the trajectory
750  NekDouble zp, velmag0, velmag, length0, zcutoff=0.0, zmax;
751 
752  zp = fabs(x2[j]);
753 
754  zmax = Vmath::Vmax(nq, x2, 1);
755  velmag = ComputeCirculatingArclength(x2[j], radius*radius);
756 
757  switch(m_surfaceType)
758  {
759 
760  case SolverUtils::eSphere:
761  case SolverUtils::eTRSphere:
762  {
763  zcutoff = 0.7;
764  velmag0 = velmag;
765  }
766  break;
767 
768  case SolverUtils::eIrregular:
769  {
770  zcutoff = 0.7;
771  length0 = 0.88;
772  // velmag0 = length0*( 1.0 - (zp -
773 zcutoff)/(1.0-zcutoff) );
774  velmag0 = (length0/(zcutoff*zcutoff-zmax*zmax))*( zp*zp -
775 zmax*zmax);
776  }
777  break;
778 
779  case SolverUtils::eNonconvex:
780  {
781  zcutoff = 0.7;
782  length0 = 1.21;
783  // velmag0 = length0*( 1.0 - (zp -
784 zcutoff)/(1.0-zcutoff) );
785  velmag0 = (length0/(zcutoff*zcutoff - 1.0))*( zp*zp -
786 1.0);
787  }
788  break;
789 
790  default:
791  break;
792  }
793 
794  if( zp > zcutoff )
795  {
796  velmag = velmag0;
797  }
798 
799  vel_phi = m_waveFreq*velmag;
800  vel_theta = 0.0;
801 
802  xyrad = sqrt( x0[j]*x0[j] + x1[j]*x1[j] );
803  if(xyrad<Tol)
804  {
805  sin_phi = 0.0;
806  cos_phi = 0.0;
807  }
808 
809  else
810  {
811  sin_phi = x1[j]/xyrad;
812  cos_phi = x0[j]/xyrad;
813  }
814 
815  // sin_theta = x2[j]/radius;
816  cos_theta = sqrt( x0[j]*x0[j] + x1[j]*x1[j] )/radius;
817 
818  if(fabs(velmag) < 0.001)
819  {
820  velocity[0][j] = 0.0 ;
821  velocity[1][j] = 0.0 ;
822  velocity[2][j] = 0.0;
823  }
824 
825  else
826  {
827  velocity[0][j] = -1.0*vel_phi*sin_phi ;
828  velocity[1][j] = vel_phi*cos_phi ;
829  velocity[2][j] = 0.0;
830  }
831 
832  } // else-loop
833  }
834 
835  // Project the veloicty on the tangent plane
836  int nq = m_fields[0]->GetNpoints();
837 
838  // Check MovingFrames \cdot SurfaceNormals = 0
839  Array<OneD, NekDouble> temp0(nq,0.0);
840  Array<OneD, NekDouble> temp1(nq,0.0);
841  Array<OneD, Array<OneD, NekDouble> > newvelocity(m_spacedim);
842 
843  for(int k=0; k<m_spacedim; ++k)
844  {
845  newvelocity[k] = Array<OneD, NekDouble>(nq);
846  }
847 
848  std::cout <<"=====================================================" <<
849 std::endl;
850  std::cout << "Velocity vector is projected onto the tangent plane " <<
851 std::endl;
852  std::cout <<"=====================================================" <<
853 std::endl;
854  GramSchumitz(m_surfaceNormal, m_velocity, newvelocity);
855 
856  for(int k=0; k<m_spacedim; ++k)
857  {
858  Vmath::Vcopy(nq, &newvelocity[k][0], 1, &m_velocity[k][0], 1);
859  }
860  }
861 
862 */
863 
865  const NekDouble Rhs)
866 {
867 
868  NekDouble Tol = 0.0001, Maxiter = 1000, N = 100;
869  NekDouble newy, F = 0.0, dF = 1.0, y0, tmp;
870 
871  Array<OneD, NekDouble> xp(N + 1);
872  Array<OneD, NekDouble> yp(N + 1);
873 
874  NekDouble intval = 0.0;
875  switch (m_surfaceType)
876  {
879  {
880  intval = sqrt(Rhs - zlevel * zlevel);
881  }
882  break;
883 
885  {
886  intval = sqrt(0.5 * (Rhs - zlevel * zlevel * zlevel * zlevel -
887  zlevel * zlevel));
888  }
889  break;
890 
892  {
893  tmp = 0.5 *
894  (Rhs - zlevel * zlevel * zlevel * zlevel - zlevel * zlevel);
895  intval = sqrt(0.5 * (1.0 + sqrt(1.0 + 4.0 * tmp)));
896  }
897  break;
898 
899  default:
900  break;
901  }
902 
903  switch (m_surfaceType)
904  {
905  // Find the half of all the xp and yp on zlevel ....
909  {
910  for (int j = 0; j < N + 1; ++j)
911  {
912  xp[j] = j * 2.0 * intval / N - intval;
913 
914  y0 = 1.0;
915  for (int i = 0; i < Maxiter; ++i)
916  {
917  switch (m_surfaceType)
918  {
919  // Find the half of all the xp and yp on zlevel ....
922  {
923  F = xp[j] * xp[j] + y0 * y0 + zlevel * zlevel - Rhs;
924  dF = 2.0 * y0;
925  }
926  break;
927 
929  {
930  F = 2.0 * xp[j] * xp[j] + y0 * y0 * y0 * y0 +
931  y0 * y0 + zlevel * zlevel * zlevel * zlevel +
932  zlevel * zlevel - Rhs;
933  dF = 4.0 * y0 * y0 * y0 + 2.0 * y0;
934  }
935  break;
936 
937  default:
938  break;
939  }
940 
941  newy = y0 - F / dF;
942 
943  if (fabs(F / dF) < Tol)
944  {
945  yp[j] = newy;
946  break;
947  }
948 
949  else
950  {
951  y0 = newy;
952  }
953 
954  ASSERTL0(i < Maxiter,
955  "Advection Velocity convergence fails");
956 
957  } // i-loop
958  }
959  }
960  break;
961 
963  {
964  for (int j = 0; j < N + 1; ++j)
965  {
966  xp[j] = j * 2.0 * intval / N - intval;
967  tmp = 0.5 * Rhs -
968  0.5 * (zlevel * zlevel * zlevel * zlevel +
969  zlevel * zlevel) -
970  (xp[j] * xp[j] * xp[j] * xp[j] - xp[j] * xp[j]);
971  if (tmp < 0)
972  {
973  tmp = -1.0 * tmp;
974  }
975  yp[j] = sqrt(tmp);
976  } // j-loop
977  }
978  break;
979 
980  default:
981  break;
982 
983  } // switch-loop
984 
985  NekDouble pi = 3.14159265358979323846;
986  NekDouble arclength = 0.0;
987  for (int j = 0; j < N; ++j)
988  {
989  arclength =
990  arclength + sqrt((yp[j + 1] - yp[j]) * (yp[j + 1] - yp[j]) +
991  (xp[j + 1] - xp[j]) * (xp[j + 1] - xp[j])) /
992  pi;
993  }
994 
995  return arclength;
996 }
997 
999  bool dumpInitialConditions,
1000  const int domain)
1001 {
1002  boost::ignore_unused(domain);
1003 
1004  int nq = m_fields[0]->GetNpoints();
1005 
1006  Array<OneD, NekDouble> u(nq);
1007 
1008  switch (m_TestType)
1009  {
1010  case eAdvectionBell:
1011  {
1013  m_fields[0]->SetPhys(u);
1014 
1015  m_Mass0 = m_fields[0]->Integral(u);
1016 
1017  // forward transform to fill the modal coeffs
1018  for (int i = 0; i < m_fields.size(); ++i)
1019  {
1020  m_fields[i]->SetPhysState(true);
1021  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1022  m_fields[i]->UpdateCoeffs());
1023  }
1024  }
1025  break;
1026 
1027  case eTestPlane:
1028  {
1029  Test2Dproblem(initialtime, u);
1030  m_fields[0]->SetPhys(u);
1031 
1032  m_Mass0 = m_fields[0]->Integral(u);
1033 
1034  // forward transform to fill the modal coeffs
1035  for (int i = 0; i < m_fields.size(); ++i)
1036  {
1037  m_fields[i]->SetPhysState(true);
1038  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1039  m_fields[i]->UpdateCoeffs());
1040  }
1041  }
1042  break;
1043 
1044  case eTestPlaneMassConsv:
1045  {
1046  AdvectionBellPlane(u);
1047  m_fields[0]->SetPhys(u);
1048 
1049  m_Mass0 = m_fields[0]->Integral(u);
1050  std::cout << "m_Mass0 = " << m_Mass0 << std::endl;
1051 
1052  // forward transform to fill the modal coeffs
1053  for (int i = 0; i < m_fields.size(); ++i)
1054  {
1055  m_fields[i]->SetPhysState(true);
1056  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1057  m_fields[i]->UpdateCoeffs());
1058  }
1059  }
1060  break;
1061 
1062  case eTestCube:
1063  {
1064  Test3Dproblem(initialtime, u);
1065  m_fields[0]->SetPhys(u);
1066 
1067  // forward transform to fill the modal coeffs
1068  for (int i = 0; i < m_fields.size(); ++i)
1069  {
1070  m_fields[i]->SetPhysState(true);
1071  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1072  m_fields[i]->UpdateCoeffs());
1073  }
1074  }
1075  break;
1076 
1077  default:
1078  break;
1079  }
1080 
1081  if (dumpInitialConditions)
1082  {
1083  // dump initial conditions to file
1084  std::string outname = m_sessionName + "_initial.chk";
1085  WriteFld(outname);
1086  }
1087 }
1088 
1090 {
1091  int nq = m_fields[0]->GetNpoints();
1092 
1093  NekDouble dist, m_radius_limit;
1094 
1095  Array<OneD, NekDouble> x(nq);
1096  Array<OneD, NekDouble> y(nq);
1097  Array<OneD, NekDouble> z(nq);
1098 
1099  m_fields[0]->GetCoords(x, y, z);
1100 
1101  // Sets of parameters
1102  m_radius_limit = 0.5;
1103 
1104  NekDouble x0j, x1j;
1105  outfield = Array<OneD, NekDouble>(nq, 0.0);
1106  for (int j = 0; j < nq; ++j)
1107  {
1108  x0j = x[j];
1109  x1j = y[j];
1110 
1111  dist = sqrt(x0j * x0j + x1j * x1j);
1112 
1113  if (dist < m_radius_limit)
1114  {
1115  outfield[j] = 0.5 * (1.0 + cos(m_pi * dist / m_radius_limit));
1116  }
1117  else
1118  {
1119  outfield[j] = 0.0;
1120  }
1121  }
1122 }
1123 
1125 {
1126  int nq = m_fields[0]->GetNpoints();
1127 
1128  NekDouble dist, radius, cosdiff, sin_theta, cos_theta, sin_varphi,
1129  cos_varphi;
1130  NekDouble m_theta_c, m_varphi_c, m_radius_limit, m_c0;
1131 
1132  Array<OneD, NekDouble> x(nq);
1133  Array<OneD, NekDouble> y(nq);
1134  Array<OneD, NekDouble> z(nq);
1135 
1136  m_fields[0]->GetCoords(x, y, z);
1137 
1138  // Sets of parameters
1139  m_theta_c = 0.0;
1140  m_varphi_c = 3.0 * m_pi / 2.0;
1141  m_radius_limit = 7.0 * m_pi / 64.0;
1142  m_c0 = 0.0;
1143 
1144  NekDouble x0j, x1j, x2j;
1145  outfield = Array<OneD, NekDouble>(nq, 0.0);
1146  for (int j = 0; j < nq; ++j)
1147  {
1148  x0j = x[j];
1149  x1j = y[j];
1150  x2j = z[j];
1151 
1152  radius = sqrt(x0j * x0j + x1j * x1j + x2j * x2j);
1153 
1154  sin_varphi = x1j / sqrt(x0j * x0j + x1j * x1j);
1155  cos_varphi = x0j / sqrt(x0j * x0j + x1j * x1j);
1156 
1157  sin_theta = x2j / radius;
1158  cos_theta = sqrt(x0j * x0j + x1j * x1j) / radius;
1159 
1160  cosdiff = cos_varphi * cos(m_varphi_c) + sin_varphi * sin(m_varphi_c);
1161  dist = radius * acos(sin(m_theta_c) * sin_theta +
1162  cos(m_theta_c) * cos_theta * cosdiff);
1163 
1164  if (dist < m_radius_limit)
1165  {
1166  outfield[j] =
1167  0.5 * (1.0 + cos(m_pi * dist / m_radius_limit)) + m_c0;
1168  }
1169  else
1170  {
1171  outfield[j] = m_c0;
1172  }
1173  }
1174 }
1175 
1177  Array<OneD, NekDouble> &outfield)
1178 {
1179  int nq = m_fields[0]->GetNpoints();
1180 
1181  Array<OneD, NekDouble> x0(nq);
1182  Array<OneD, NekDouble> x1(nq);
1183  Array<OneD, NekDouble> x2(nq);
1184 
1185  m_fields[0]->GetCoords(x0, x1, x2);
1186 
1187  Array<OneD, NekDouble> u(nq);
1188 
1189  for (int i = 0; i < nq; ++i)
1190  {
1191  u[i] = cos(m_pi * (x0[i] - m_advx * time)) *
1192  cos(m_pi * (x1[i] - m_advy * time));
1193  }
1194 
1195  outfield = u;
1196 }
1197 
1199  Array<OneD, NekDouble> &outfield)
1200 {
1201  int nq = m_fields[0]->GetNpoints();
1202 
1203  Array<OneD, NekDouble> x0(nq);
1204  Array<OneD, NekDouble> x1(nq);
1205  Array<OneD, NekDouble> x2(nq);
1206 
1207  m_fields[0]->GetCoords(x0, x1, x2);
1208 
1209  Array<OneD, NekDouble> u(nq);
1210 
1211  for (int i = 0; i < nq; ++i)
1212  {
1213  u[i] = cos(m_pi * (x0[i] - m_advx * time)) *
1214  cos(m_pi * (x1[i] - m_advy * time)) *
1215  cos(m_pi * (x2[i] - m_advz * time));
1216  }
1217 
1218  outfield = u;
1219 }
1220 
1222 {
1223  int nq = m_fields[0]->GetNpoints();
1224 
1225  Array<OneD, NekDouble> velcoeff(nq, 0.0);
1226 
1227  Array<OneD, NekDouble> Dtmp0(nq);
1228  Array<OneD, NekDouble> Dtmp1(nq);
1229  Array<OneD, NekDouble> Dtmp2(nq);
1230  Array<OneD, NekDouble> Drv(nq);
1231 
1232  vellc = Array<OneD, NekDouble>(nq, 0.0);
1233 
1234  // m_vellc = \nabla m_vel \cdot tan_i
1235  Array<OneD, NekDouble> tmp(nq);
1236  Array<OneD, NekDouble> vessel(nq);
1237 
1238  for (int j = 0; j < m_shapedim; ++j)
1239  {
1240  Vmath::Zero(nq, velcoeff, 1);
1241  for (int k = 0; k < m_spacedim; ++k)
1242  {
1243  // a_j = tan_j cdot m_vel
1244  Vmath::Vvtvp(nq, &m_movingframes[j][k * nq], 1, &m_velocity[k][0],
1245  1, &velcoeff[0], 1, &velcoeff[0], 1);
1246  }
1247 
1248  // d a_j / d x^k
1249  m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
1250 
1251  for (int k = 0; k < m_spacedim; ++k)
1252  {
1253  // tan_j^k ( d a_j / d x^k )
1254  switch (k)
1255  {
1256  case (0):
1257  {
1258  Vmath::Vvtvp(nq, &Dtmp0[0], 1, &m_movingframes[j][k * nq],
1259  1, &vellc[0], 1, &vellc[0], 1);
1260  }
1261  break;
1262 
1263  case (1):
1264  {
1265  Vmath::Vvtvp(nq, &Dtmp1[0], 1, &m_movingframes[j][k * nq],
1266  1, &vellc[0], 1, &vellc[0], 1);
1267  }
1268  break;
1269 
1270  case (2):
1271  {
1272  Vmath::Vvtvp(nq, &Dtmp2[0], 1, &m_movingframes[j][k * nq],
1273  1, &vellc[0], 1, &vellc[0], 1);
1274  }
1275  break;
1276  }
1277  }
1278  }
1279 }
1280 
1282  Array<OneD, Array<OneD, NekDouble>> &veldotMF)
1283 {
1284  int nq = m_fields[0]->GetNpoints();
1285 
1287 
1288  Array<OneD, NekDouble> magMF(nq);
1289  for (int j = 0; j < m_shapedim; ++j)
1290  {
1291  veldotMF[j] = Array<OneD, NekDouble>(nq, 0.0);
1292 
1293  for (int k = 0; k < m_spacedim; ++k)
1294  {
1295  Vmath::Vvtvp(nq, &m_movingframes[j][k * nq], 1, &m_velocity[k][0],
1296  1, &veldotMF[j][0], 1, &veldotMF[j][0], 1);
1297  }
1298  }
1299 }
1300 
1302  Array<OneD, NekDouble> &outfield,
1303  const NekDouble time)
1304 {
1305  boost::ignore_unused(field);
1306 
1307  switch (m_TestType)
1308  {
1309  case eAdvectionBell:
1310  {
1311  AdvectionBellSphere(outfield);
1312  }
1313  break;
1314 
1315  case eTestPlane:
1316  {
1317  Test2Dproblem(time, outfield);
1318  }
1319  break;
1320 
1321  case eTestPlaneMassConsv:
1322  {
1323  AdvectionBellPlane(outfield);
1324  }
1325  break;
1326 
1327  case eTestCube:
1328  {
1329  Test3Dproblem(time, outfield);
1330  }
1331  break;
1332 
1333  default:
1334  break;
1335  }
1336 }
1337 
1339 {
1342  SolverUtils::AddSummaryItem(s, "Rotation Angle", m_RotAngle);
1343 }
1344 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
@ eTestPlaneMassConsv
Definition: MMFAdvection.h:46
@ eAdvectionBell
Definition: MMFAdvection.h:48
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
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:68
void Test2Dproblem(const NekDouble time, Array< OneD, NekDouble > &outfield)
virtual void v_DoSolve()
Solves an unsteady problem.
virtual ~MMFAdvection()
Destructor.
Array< OneD, Array< OneD, NekDouble > > m_veldotMF
Definition: MMFAdvection.h:98
void Test3Dproblem(const NekDouble time, Array< OneD, NekDouble > &outfield)
void AdvectionBellSphere(Array< OneD, NekDouble > &outfield)
NekDouble ComputeCirculatingArclength(const NekDouble zlevel, const NekDouble Rhs)
MMFAdvection(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Session reader.
void ComputeNablaCdotVelocity(Array< OneD, NekDouble > &vellc)
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print Summary.
void EvaluateAdvectionVelocity(Array< OneD, Array< OneD, NekDouble >> &velocity)
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the RHS.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the projection.
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
void WeakDGDirectionalAdvection(const Array< OneD, Array< OneD, NekDouble >> &InField, Array< OneD, Array< OneD, NekDouble >> &OutField)
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble >>> &flux)
Evaluate the flux at each solution point.
void ComputeveldotMF(Array< OneD, Array< OneD, NekDouble >> &veldotMF)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFAdvection.h:68
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Definition: MMFAdvection.h:86
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
Definition: MMFAdvection.h:95
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
virtual void v_SetInitialConditions(const NekDouble initialtime, bool dumpInitialConditions, const int domain)
static std::string className
Name of class.
Definition: MMFAdvection.h:78
Array< OneD, NekDouble > m_vellc
Definition: MMFAdvection.h:99
void AdvectionBellPlane(Array< OneD, NekDouble > &outfield)
virtual void v_InitObject(bool DeclareFields=true)
Initialise the object.
Array< OneD, NekDouble > m_traceVn
Definition: MMFAdvection.h:96
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
int m_spacedim
Spatial dimension (>= expansion dim).
SOLVER_UTILS_EXPORT int GetTraceNpoints()
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
NekDouble m_fintime
Finish time of the simulation.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
NekDouble m_checktime
Time between checkpoints.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
SOLVER_UTILS_EXPORT int GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
int m_steps
Number of steps to take.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
A base class for PDEs which include an advection component.
Definition: MMFSystem.h:147
SOLVER_UTILS_EXPORT void GramSchumitz(const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &outarray, bool KeepTheMagnitude=true)
Definition: MMFSystem.cpp:2405
SOLVER_UTILS_EXPORT void VectorCrossProd(const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &v3)
Definition: MMFSystem.cpp:699
SOLVER_UTILS_EXPORT NekDouble AvgInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2298
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:186
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:53
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:190
SOLVER_UTILS_EXPORT void ComputencdotMF()
Definition: MMFSystem.cpp:320
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2344
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:189
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2469
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
int m_infosteps
Number of time steps between outputting status information.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekZeroTol
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:47
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:48
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
@ eTestPlane
Definition: MMFDiffusion.h:48
@ eTestCube
Definition: MMFDiffusion.h:49
@ SIZE_TestType
Length of enum list.
Definition: MMFDiffusion.h:55
const char *const TestTypeMap[]
Definition: MMFDiffusion.h:58
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 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
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291