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
43
44namespace Nektar
45{
46
47std::string MMFAdvection::className =
49 "MMFAdvection", MMFAdvection::create, "MMFAdvection equation.");
50
53 : UnsteadySystem(pSession, pGraph), MMFSystem(pSession, pGraph),
54 AdvectionSystem(pSession, pGraph)
55{
56}
57
58/**
59 * @brief Initialisation object for the unsteady linear advection equation.
60 */
61void MMFAdvection::v_InitObject(bool DeclareFields)
62{
63 // Call to the initialisation object
64 UnsteadySystem::v_InitObject(DeclareFields);
65
66 int nq = m_fields[0]->GetNpoints();
67 int shapedim = m_fields[0]->GetShapeDimension();
68 Array<OneD, Array<OneD, NekDouble>> Anisotropy(shapedim);
69 for (int j = 0; j < shapedim; ++j)
70 {
71 Anisotropy[j] = Array<OneD, NekDouble>(nq, 1.0);
72 }
73
74 MMFSystem::MMFInitObject(Anisotropy);
75
76 // Define TestType
77 ASSERTL0(m_session->DefinesSolverInfo("TESTTYPE"),
78 "No TESTTYPE defined in session.");
79 std::string TestTypeStr = m_session->GetSolverInfo("TESTTYPE");
80 for (int i = 0; i < (int)SIZE_TestType; ++i)
81 {
82 if (boost::iequals(TestTypeMap[i], TestTypeStr))
83 {
85 break;
86 }
87 }
88
89 m_session->LoadParameter("Angular Frequency", m_waveFreq, m_pi);
90 m_session->LoadParameter("Rotational Angle", m_RotAngle, 0.0);
91 m_session->LoadParameter("Velocity Projection", m_VelProjection, 0);
92
93 // Read the advection velocities from session file
94 m_session->LoadParameter("advx", m_advx, 1.0);
95 m_session->LoadParameter("advy", m_advy, 1.0);
96 m_session->LoadParameter("advz", m_advz, 1.0);
97
98 std::vector<std::string> vel;
99 vel.push_back("Vx");
100 vel.push_back("Vy");
101 vel.push_back("Vz");
102
103 // Resize the advection velocities vector to dimension of the problem
104 vel.resize(m_spacedim);
105
106 // Store in the global variable m_velocity the advection velocities
108 for (int k = 0; k < m_spacedim; ++k)
109 {
111 }
112
113 switch (m_surfaceType)
114 {
119 {
120 // true = project velocity onto the tangent plane
122 }
123 break;
124
127 {
128 GetFunction("AdvectionVelocity")->Evaluate(vel, m_velocity);
129 }
130 break;
131
132 default:
133 break;
134 }
135
136 std::cout << "|Velocity vector| = ( " << RootMeanSquare(m_velocity[0])
137 << " , " << RootMeanSquare(m_velocity[1]) << " , "
138 << RootMeanSquare(m_velocity[2]) << " ) " << std::endl;
139
140 // Define the normal velocity fields
141 if (m_fields[0]->GetTrace())
142 {
144 }
145
146 std::string advName;
147 std::string riemName;
148 m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
151 m_advObject->SetFluxVector(&MMFAdvection::GetFluxVector, this);
152 m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
154 riemName, m_session);
155 m_riemannSolver->SetScalar("Vn", &MMFAdvection::GetNormalVelocity, this);
156
157 m_advObject->SetRiemannSolver(m_riemannSolver);
158 m_advObject->InitObject(m_session, m_fields);
159
160 // Compute m_traceVn = n \cdot v
162
163 // Compute m_vellc = nabal a^j \cdot m_vel
165 std::cout << "m_vellc is generated with mag = " << AvgInt(m_vellc)
166 << std::endl;
167
168 // Compute vel \cdot MF
170
171 // Modify e^i as v^i e^i
172 for (int j = 0; j < m_shapedim; j++)
173 {
174 for (int k = 0; k < m_spacedim; k++)
175 {
176 Vmath::Vmul(nq, &m_veldotMF[j][0], 1, &m_movingframes[j][k * nq], 1,
177 &m_movingframes[j][k * nq], 1);
178 }
179 }
180
181 // Reflect it into m_ncdotMFFwd and Bwd
183
184 // If explicit it computes RHS and PROJECTION for the time integration
186 {
189 }
190 // Otherwise it gives an error (no implicit integration)
191 else
192 {
193 ASSERTL0(false, "Implicit unsteady Advection not set up.");
194 }
195}
196
198{
199 ASSERTL0(m_intScheme != nullptr, "No time integration scheme.");
200
201 int i, nchk = 1;
202 int nvariables = 0;
203 int nfields = m_fields.size();
204 int nq = m_fields[0]->GetNpoints();
205
206 if (m_intVariables.empty())
207 {
208 for (i = 0; i < nfields; ++i)
209 {
210 m_intVariables.push_back(i);
211 }
212 nvariables = nfields;
213 }
214 else
215 {
216 nvariables = m_intVariables.size();
217 }
218
219 // Set up wrapper to fields data storage.
220 Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
221 Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
222
223 // Order storage to list time-integrated fields first.
224 for (i = 0; i < nvariables; ++i)
225 {
226 fields[i] = m_fields[m_intVariables[i]]->GetPhys();
227 m_fields[m_intVariables[i]]->SetPhysState(false);
228 }
229
230 // Initialise time integration scheme
231 m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
232
233 // Check uniqueness of checkpoint output
234 ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
235 (m_checktime > 0.0 && m_checksteps == 0) ||
236 (m_checktime == 0.0 && m_checksteps > 0),
237 "Only one of IO_CheckTime and IO_CheckSteps "
238 "should be set!");
239
241 bool doCheckTime = false;
242 int step = 0;
243 NekDouble intTime = 0.0;
244 NekDouble cpuTime = 0.0;
245 NekDouble elapsed = 0.0;
246
247 int Ntot, indx;
248 // Perform integration in time.
249 Ntot = m_checksteps ? m_steps / m_checksteps + 1 : 0;
250
251 Array<OneD, NekDouble> dMass(Ntot ? Ntot : 1);
252
253 Array<OneD, NekDouble> zeta(nq);
254 Array<OneD, Array<OneD, NekDouble>> fieldsprimitive(nvariables);
255 for (int i = 0; i < nvariables; ++i)
256 {
257 fieldsprimitive[i] = Array<OneD, NekDouble>(nq);
258 }
259
261 {
262 timer.Start();
263 fields = m_intScheme->TimeIntegrate(step, m_timestep);
264 timer.Stop();
265
267 elapsed = timer.TimePerTest(1);
268 intTime += elapsed;
269 cpuTime += elapsed;
270
271 // Write out status information
272 if (m_infosteps && !((step + 1) % m_infosteps) &&
273 m_session->GetComm()->GetRank() == 0)
274 {
275 std::cout << "Steps: " << std::setw(8) << std::left << step + 1
276 << " "
277 << "Time: " << std::setw(12) << std::left << m_time;
278
279 std::stringstream ss;
280 ss << cpuTime << "s";
281 std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
282 << std::endl;
283
284 // Masss = h^*
285 indx = m_checksteps ? (step + 1) / m_checksteps : 0;
286 dMass[indx] =
287 (m_fields[0]->Integral(fields[0]) - m_Mass0) / m_Mass0;
288
289 std::cout << "dMass = " << std::setw(8) << std::left << dMass[indx]
290 << std::endl;
291
292 cpuTime = 0.0;
293 }
294
295 // Transform data into coefficient space
296 for (i = 0; i < nvariables; ++i)
297 {
298 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
299 m_fields[m_intVariables[i]]->FwdTransLocalElmt(
300 fields[i], m_fields[m_intVariables[i]]->UpdateCoeffs());
301 m_fields[m_intVariables[i]]->SetPhysState(false);
302 }
303
304 // Write out checkpoint files
305 if ((m_checksteps && step && !((step + 1) % m_checksteps)) ||
306 doCheckTime)
307 {
308 Checkpoint_Output(nchk++);
309 doCheckTime = false;
310 }
311
312 // Step advance
313 ++step;
314 }
315
316 std::cout << "dMass = ";
317 for (i = 0; i < Ntot; ++i)
318 {
319 std::cout << dMass[i] << " , ";
320 }
321 std::cout << std::endl << std::endl;
322
323 // Print out summary statistics
324 if (m_session->GetComm()->GetRank() == 0)
325 {
326 if (m_cflSafetyFactor > 0.0)
327 {
328 std::cout << "CFL safety factor : " << m_cflSafetyFactor
329 << std::endl
330 << "CFL time-step : " << m_timestep << std::endl;
331 }
332
333 if (m_session->GetSolverInfo("Driver") != "SteadyState")
334 {
335 std::cout << "Time-integration : " << intTime << "s" << std::endl;
336 }
337 }
338
339 for (i = 0; i < nvariables; ++i)
340 {
341 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
342 m_fields[m_intVariables[i]]->SetPhysState(true);
343 }
344
345 for (i = 0; i < nvariables; ++i)
346 {
347 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
348 m_fields[i]->UpdateCoeffs());
349 }
350}
351
352/**
353 * @brief Get the normal velocity for the linear advection equation.
354 */
356{
357 // Number of trace (interface) points
358 int i;
359 int nTracePts = GetTraceNpoints();
360
361 // Auxiliary variable to compute the normal velocity
362 Array<OneD, NekDouble> tmp(nTracePts);
363
364 // Reset the normal velocity
365 Vmath::Zero(nTracePts, m_traceVn, 1);
366
367 for (i = 0; i < m_velocity.size(); ++i)
368 {
369 m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
370
371 Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, tmp, 1, m_traceVn, 1,
372 m_traceVn, 1);
373 }
374
375 return m_traceVn;
376}
377
378/**
379 * @brief Compute the right-hand side for the linear advection equation.
380 *
381 * @param inarray Given fields.
382 * @param outarray Calculated solution.
383 * @param time Time.
384 */
386 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
388 [[maybe_unused]] const NekDouble time)
389{
390 int i;
391 int nvariables = inarray.size();
392 int npoints = GetNpoints();
393
394 switch (m_projectionType)
395 {
397 {
398 int ncoeffs = inarray[0].size();
399
400 if (m_spacedim == 3)
401 {
402 Array<OneD, Array<OneD, NekDouble>> WeakAdv(nvariables);
403
404 WeakAdv[0] = Array<OneD, NekDouble>(ncoeffs * nvariables);
405 for (i = 1; i < nvariables; ++i)
406 {
407 WeakAdv[i] = WeakAdv[i - 1] + ncoeffs;
408 }
409
410 // Compute \nabla \cdot \vel u according to MMF scheme
411 WeakDGDirectionalAdvection(inarray, WeakAdv);
412
413 for (i = 0; i < nvariables; ++i)
414 {
415 m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
416 m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
417 Vmath::Neg(npoints, outarray[i], 1);
418 }
419 }
420 else
421 {
422 m_advObject->Advect(2, m_fields, m_velocity, inarray, outarray,
423 0.0);
424
425 for (i = 0; i < nvariables; ++i)
426 {
427 Vmath::Neg(npoints, outarray[i], 1);
428 }
429 }
430 }
431 break;
432
433 default:
434 {
435 ASSERTL0(false, "Unknown projection scheme");
436 }
437 break;
438 }
439}
440
441/**
442 * @brief Compute the projection for the linear advection equation.
443 *
444 * @param inarray Given fields.
445 * @param outarray Calculated solution.
446 * @param time Time.
447 */
449 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
450 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
451{
452 // Counter variable
453 int i;
454
455 // Number of fields (variables of the problem)
456 int nVariables = inarray.size();
457
458 // Set the boundary conditions
460
461 // Switch on the projection type (Discontinuous or Continuous)
462 switch (m_projectionType)
463 {
464 // Discontinuous projection
466 {
467 // Number of quadrature points
468 int nQuadraturePts = GetNpoints();
469
470 // Just copy over array
471 if (inarray != outarray)
472 {
473 for (i = 0; i < nVariables; ++i)
474 {
475 Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
476 }
477 }
478 break;
479 }
480
481 // Continuous projection
484 {
485 Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs(), 0.0);
486 for (i = 0; i < nVariables; ++i)
487 {
488 m_fields[i]->FwdTrans(inarray[i], coeffs);
489 m_fields[i]->BwdTrans(coeffs, outarray[i]);
490 }
491 break;
492 }
493
494 default:
495 ASSERTL0(false, "Unknown projection scheme");
496 break;
497 }
498}
499
500/**
501 * @brief Return the flux vector for the linear advection equation.
502 *
503 * @param i Component of the flux vector to calculate.
504 * @param physfield Fields.
505 * @param flux Resulting flux.
506 */
508 const Array<OneD, Array<OneD, NekDouble>> &physfield,
510{
511 ASSERTL1(flux[0].size() == m_velocity.size(),
512 "Dimension of flux array and velocity array do not match");
513
514 int i, j;
515 int nq = physfield[0].size();
516
517 for (i = 0; i < flux.size(); ++i)
518 {
519 for (j = 0; j < flux[0].size(); ++j)
520 {
521 Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1, flux[i][j], 1);
522 }
523 }
524}
525
527 const Array<OneD, Array<OneD, NekDouble>> &InField,
529{
530 int i, j;
531 int ncoeffs = GetNcoeffs();
532 int nTracePointsTot = GetTraceNpoints();
533 int nvariables = m_fields.size();
534
535 Array<OneD, Array<OneD, NekDouble>> physfield(nvariables);
536
537 // Get the variables in physical space
538 // already in physical space
539 for (i = 0; i < nvariables; ++i)
540 {
541 physfield[i] = InField[i];
542 }
543
545 for (i = 0; i < nvariables; ++i)
546 {
547 for (j = 0; j < m_shapedim; ++j)
548 {
549 WeakDeriv[j] = Array<OneD, NekDouble>(ncoeffs, 0.0);
550
551 // Directional derivation with respect to the j'th moving frame
552 // tmp[j] = \nabla \physfield[i] \cdot \mathbf{e}^j
553 // Implemented at TriExp::v_IProductWRTDirectionalDerivBase_SumFa
554 m_fields[0]->IProductWRTDirectionalDerivBase(
555 m_movingframes[j], physfield[i], WeakDeriv[j]);
556 }
557
558 // Get the numerical flux and add to the modal coeffs
559 // if the NumericalFluxs function already includes the
560 // normal in the output
561
562 Array<OneD, NekDouble> Fwd(nTracePointsTot);
563 Array<OneD, NekDouble> Bwd(nTracePointsTot);
564
565 Array<OneD, NekDouble> flux(nTracePointsTot, 0.0);
566 Array<OneD, NekDouble> fluxFwd(nTracePointsTot);
567 Array<OneD, NekDouble> fluxBwd(nTracePointsTot);
568
569 // Evaluate numerical flux in physical space which may in
570 // general couple all component of vectors
571 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
572
573 // evaulate upwinded m_fields[i]
574 m_fields[i]->GetTrace()->Upwind(m_traceVn, Fwd, Bwd, flux);
575
576 OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
577 for (j = 0; j < m_shapedim; ++j)
578 {
579 // calculate numflux = (n \cdot MF)*flux
580 Vmath::Vmul(nTracePointsTot, &flux[0], 1, &m_ncdotMFFwd[j][0], 1,
581 &fluxFwd[0], 1);
582 Vmath::Vmul(nTracePointsTot, &flux[0], 1, &m_ncdotMFBwd[j][0], 1,
583 &fluxBwd[0], 1);
584 Vmath::Neg(ncoeffs, WeakDeriv[j], 1);
585
586 // FwdBwdtegral because generallize (N \cdot MF)_{FWD} \neq -(N
587 // \cdot MF)_{BWD}
588 m_fields[i]->AddFwdBwdTraceIntegral(fluxFwd, fluxBwd, WeakDeriv[j]);
589 m_fields[i]->SetPhysState(false);
590
591 Vmath::Vadd(ncoeffs, &WeakDeriv[j][0], 1, &OutField[i][0], 1,
592 &OutField[i][0], 1);
593 }
594 }
595}
596
599{
600 int nq = m_fields[0]->GetNpoints();
601
602 NekDouble vel_phi, vel_theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
603 NekDouble x0j, x1j, x2j;
604
608
609 m_fields[0]->GetCoords(x0, x1, x2);
610
611 // theta = a*sin(z/r), phi = a*tan(y/x);
612 for (int j = 0; j < nq; j++)
613 {
614 x0j = x0[j];
615 x1j = x1[j];
616 x2j = x2[j];
617
618 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
619 cos_theta);
620
621 vel_phi = m_waveFreq * (cos_theta * cos(m_RotAngle) +
622 sin_theta * cos_varphi * sin(m_RotAngle));
623 vel_theta = -1.0 * m_waveFreq * sin_theta * sin(m_RotAngle);
624
625 velocity[0][j] =
626 -1.0 * vel_phi * sin_varphi - vel_theta * sin_theta * cos_varphi;
627 velocity[1][j] =
628 vel_phi * cos_varphi - vel_theta * sin_theta * sin_varphi;
629 velocity[2][j] = vel_theta * cos_theta;
630 }
631
632 // Project the veloicty on the tangent plane
633
634 if (m_VelProjection)
635 {
636 // Check MovingFrames \cdot SurfaceNormals = 0
638
642
643 for (int k = 0; k < m_spacedim; ++k)
644 {
645 newvelocity[k] = Array<OneD, NekDouble>(nq);
646 MF1[k] = Array<OneD, NekDouble>(nq);
647 MF2[k] = Array<OneD, NekDouble>(nq);
648 SN[k] = Array<OneD, NekDouble>(nq);
649
650 Vmath::Vcopy(nq, &m_movingframes[0][k * nq], 1, &MF1[k][0], 1);
651 Vmath::Vcopy(nq, &m_movingframes[1][k * nq], 1, &MF2[k][0], 1);
652 }
653
654 VectorCrossProd(MF1, MF2, SN);
655 GramSchumitz(SN, m_velocity, newvelocity, true);
656
657 Array<OneD, NekDouble> tmp(nq, 0.0);
659
660 for (int k = 0; k < m_spacedim; ++k)
661 {
662 Vmath::Vsub(nq, &m_velocity[0][0], 1, &newvelocity[0][0], 1,
663 &tmp[0], 1);
664 Vmath::Vmul(nq, &tmp[0], 1, &tmp[0], 1, &tmp[0], 1);
665 Vmath::Vadd(nq, tmp, 1, Projection, 1, Projection, 1);
666 }
667
668 std::cout
669 << "Velocity vector is projected onto the tangent plane: Diff = "
670 << RootMeanSquare(Projection) << std::endl;
671
672 for (int k = 0; k < m_spacedim; ++k)
673 {
674 Vmath::Vcopy(nq, &newvelocity[k][0], 1, &m_velocity[k][0], 1);
675 }
676 }
677}
678
680 const NekDouble Rhs)
681{
682
683 NekDouble Tol = 0.0001, Maxiter = 1000, N = 100;
684 NekDouble newy, F = 0.0, dF = 1.0, y0, tmp;
685
686 Array<OneD, NekDouble> xp(N + 1);
687 Array<OneD, NekDouble> yp(N + 1);
688
689 NekDouble intval = 0.0;
690 switch (m_surfaceType)
691 {
694 {
695 intval = sqrt(Rhs - zlevel * zlevel);
696 }
697 break;
698
700 {
701 intval = sqrt(0.5 * (Rhs - zlevel * zlevel * zlevel * zlevel -
702 zlevel * zlevel));
703 }
704 break;
705
707 {
708 tmp = 0.5 *
709 (Rhs - zlevel * zlevel * zlevel * zlevel - zlevel * zlevel);
710 intval = sqrt(0.5 * (1.0 + sqrt(1.0 + 4.0 * tmp)));
711 }
712 break;
713
714 default:
715 break;
716 }
717
718 switch (m_surfaceType)
719 {
720 // Find the half of all the xp and yp on zlevel ....
724 {
725 for (int j = 0; j < N + 1; ++j)
726 {
727 xp[j] = j * 2.0 * intval / N - intval;
728
729 y0 = 1.0;
730 for (int i = 0; i < Maxiter; ++i)
731 {
732 switch (m_surfaceType)
733 {
734 // Find the half of all the xp and yp on zlevel ....
737 {
738 F = xp[j] * xp[j] + y0 * y0 + zlevel * zlevel - Rhs;
739 dF = 2.0 * y0;
740 }
741 break;
742
744 {
745 F = 2.0 * xp[j] * xp[j] + y0 * y0 * y0 * y0 +
746 y0 * y0 + zlevel * zlevel * zlevel * zlevel +
747 zlevel * zlevel - Rhs;
748 dF = 4.0 * y0 * y0 * y0 + 2.0 * y0;
749 }
750 break;
751
752 default:
753 break;
754 }
755
756 newy = y0 - F / dF;
757
758 if (fabs(F / dF) < Tol)
759 {
760 yp[j] = newy;
761 break;
762 }
763
764 else
765 {
766 y0 = newy;
767 }
768
769 ASSERTL0(i < Maxiter,
770 "Advection Velocity convergence fails");
771
772 } // i-loop
773 }
774 }
775 break;
776
778 {
779 for (int j = 0; j < N + 1; ++j)
780 {
781 xp[j] = j * 2.0 * intval / N - intval;
782 tmp = 0.5 * Rhs -
783 0.5 * (zlevel * zlevel * zlevel * zlevel +
784 zlevel * zlevel) -
785 (xp[j] * xp[j] * xp[j] * xp[j] - xp[j] * xp[j]);
786 if (tmp < 0)
787 {
788 tmp = -1.0 * tmp;
789 }
790 yp[j] = sqrt(tmp);
791 } // j-loop
792 }
793 break;
794
795 default:
796 break;
797
798 } // switch-loop
799
800 NekDouble pi = 3.14159265358979323846;
801 NekDouble arclength = 0.0;
802 for (int j = 0; j < N; ++j)
803 {
804 arclength =
805 arclength + sqrt((yp[j + 1] - yp[j]) * (yp[j + 1] - yp[j]) +
806 (xp[j + 1] - xp[j]) * (xp[j + 1] - xp[j])) /
807 pi;
808 }
809
810 return arclength;
811}
812
814 bool dumpInitialConditions,
815 [[maybe_unused]] const int domain)
816{
817 int nq = m_fields[0]->GetNpoints();
818
820
821 switch (m_TestType)
822 {
823 case eAdvectionBell:
824 {
826 m_fields[0]->SetPhys(u);
827
828 m_Mass0 = m_fields[0]->Integral(u);
829
830 // forward transform to fill the modal coeffs
831 for (int i = 0; i < m_fields.size(); ++i)
832 {
833 m_fields[i]->SetPhysState(true);
834 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
835 m_fields[i]->UpdateCoeffs());
836 }
837 }
838 break;
839
840 case eTestPlane:
841 {
842 Test2Dproblem(initialtime, u);
843 m_fields[0]->SetPhys(u);
844
845 m_Mass0 = m_fields[0]->Integral(u);
846
847 // forward transform to fill the modal coeffs
848 for (int i = 0; i < m_fields.size(); ++i)
849 {
850 m_fields[i]->SetPhysState(true);
851 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
852 m_fields[i]->UpdateCoeffs());
853 }
854 }
855 break;
856
858 {
860 m_fields[0]->SetPhys(u);
861
862 m_Mass0 = m_fields[0]->Integral(u);
863 std::cout << "m_Mass0 = " << m_Mass0 << std::endl;
864
865 // forward transform to fill the modal coeffs
866 for (int i = 0; i < m_fields.size(); ++i)
867 {
868 m_fields[i]->SetPhysState(true);
869 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
870 m_fields[i]->UpdateCoeffs());
871 }
872 }
873 break;
874
875 case eTestCube:
876 {
877 Test3Dproblem(initialtime, u);
878 m_fields[0]->SetPhys(u);
879
880 // forward transform to fill the modal coeffs
881 for (int i = 0; i < m_fields.size(); ++i)
882 {
883 m_fields[i]->SetPhysState(true);
884 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
885 m_fields[i]->UpdateCoeffs());
886 }
887 }
888 break;
889
890 default:
891 break;
892 }
893
894 if (dumpInitialConditions)
895 {
896 // dump initial conditions to file
897 std::string outname = m_sessionName + "_initial.chk";
898 WriteFld(outname);
899 }
900}
901
903{
904 int nq = m_fields[0]->GetNpoints();
905
906 NekDouble dist, m_radius_limit;
907
911
912 m_fields[0]->GetCoords(x, y, z);
913
914 // Sets of parameters
915 m_radius_limit = 0.5;
916
917 NekDouble x0j, x1j;
918 outfield = Array<OneD, NekDouble>(nq, 0.0);
919 for (int j = 0; j < nq; ++j)
920 {
921 x0j = x[j];
922 x1j = y[j];
923
924 dist = sqrt(x0j * x0j + x1j * x1j);
925
926 if (dist < m_radius_limit)
927 {
928 outfield[j] = 0.5 * (1.0 + cos(m_pi * dist / m_radius_limit));
929 }
930 else
931 {
932 outfield[j] = 0.0;
933 }
934 }
935}
936
938{
939 int nq = m_fields[0]->GetNpoints();
940
941 NekDouble dist, radius, cosdiff, sin_theta, cos_theta, sin_varphi,
942 cos_varphi;
943 NekDouble m_theta_c, m_varphi_c, m_radius_limit, m_c0;
944
948
949 m_fields[0]->GetCoords(x, y, z);
950
951 // Sets of parameters
952 m_theta_c = 0.0;
953 m_varphi_c = 3.0 * m_pi / 2.0;
954 m_radius_limit = 7.0 * m_pi / 64.0;
955 m_c0 = 0.0;
956
957 NekDouble x0j, x1j, x2j;
958 outfield = Array<OneD, NekDouble>(nq, 0.0);
959 for (int j = 0; j < nq; ++j)
960 {
961 x0j = x[j];
962 x1j = y[j];
963 x2j = z[j];
964
965 radius = sqrt(x0j * x0j + x1j * x1j + x2j * x2j);
966
967 sin_varphi = x1j / sqrt(x0j * x0j + x1j * x1j);
968 cos_varphi = x0j / sqrt(x0j * x0j + x1j * x1j);
969
970 sin_theta = x2j / radius;
971 cos_theta = sqrt(x0j * x0j + x1j * x1j) / radius;
972
973 cosdiff = cos_varphi * cos(m_varphi_c) + sin_varphi * sin(m_varphi_c);
974 dist = radius * acos(sin(m_theta_c) * sin_theta +
975 cos(m_theta_c) * cos_theta * cosdiff);
976
977 if (dist < m_radius_limit)
978 {
979 outfield[j] =
980 0.5 * (1.0 + cos(m_pi * dist / m_radius_limit)) + m_c0;
981 }
982 else
983 {
984 outfield[j] = m_c0;
985 }
986 }
987}
988
990 Array<OneD, NekDouble> &outfield)
991{
992 int nq = m_fields[0]->GetNpoints();
993
997
998 m_fields[0]->GetCoords(x0, x1, x2);
999
1001
1002 for (int i = 0; i < nq; ++i)
1003 {
1004 u[i] = cos(m_pi * (x0[i] - m_advx * time)) *
1005 cos(m_pi * (x1[i] - m_advy * time));
1006 }
1007
1008 outfield = u;
1009}
1010
1012 Array<OneD, NekDouble> &outfield)
1013{
1014 int nq = m_fields[0]->GetNpoints();
1015
1019
1020 m_fields[0]->GetCoords(x0, x1, x2);
1021
1023
1024 for (int i = 0; i < nq; ++i)
1025 {
1026 u[i] = cos(m_pi * (x0[i] - m_advx * time)) *
1027 cos(m_pi * (x1[i] - m_advy * time)) *
1028 cos(m_pi * (x2[i] - m_advz * time));
1029 }
1030
1031 outfield = u;
1032}
1033
1035{
1036 int nq = m_fields[0]->GetNpoints();
1037
1038 Array<OneD, NekDouble> velcoeff(nq, 0.0);
1039
1040 Array<OneD, NekDouble> Dtmp0(nq);
1041 Array<OneD, NekDouble> Dtmp1(nq);
1042 Array<OneD, NekDouble> Dtmp2(nq);
1043 Array<OneD, NekDouble> Drv(nq);
1044
1045 vellc = Array<OneD, NekDouble>(nq, 0.0);
1046
1047 // m_vellc = \nabla m_vel \cdot tan_i
1048 Array<OneD, NekDouble> tmp(nq);
1049 Array<OneD, NekDouble> vessel(nq);
1050
1051 for (int j = 0; j < m_shapedim; ++j)
1052 {
1053 Vmath::Zero(nq, velcoeff, 1);
1054 for (int k = 0; k < m_spacedim; ++k)
1055 {
1056 // a_j = tan_j cdot m_vel
1057 Vmath::Vvtvp(nq, &m_movingframes[j][k * nq], 1, &m_velocity[k][0],
1058 1, &velcoeff[0], 1, &velcoeff[0], 1);
1059 }
1060
1061 // d a_j / d x^k
1062 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
1063
1064 for (int k = 0; k < m_spacedim; ++k)
1065 {
1066 // tan_j^k ( d a_j / d x^k )
1067 switch (k)
1068 {
1069 case (0):
1070 {
1071 Vmath::Vvtvp(nq, &Dtmp0[0], 1, &m_movingframes[j][k * nq],
1072 1, &vellc[0], 1, &vellc[0], 1);
1073 }
1074 break;
1075
1076 case (1):
1077 {
1078 Vmath::Vvtvp(nq, &Dtmp1[0], 1, &m_movingframes[j][k * nq],
1079 1, &vellc[0], 1, &vellc[0], 1);
1080 }
1081 break;
1082
1083 case (2):
1084 {
1085 Vmath::Vvtvp(nq, &Dtmp2[0], 1, &m_movingframes[j][k * nq],
1086 1, &vellc[0], 1, &vellc[0], 1);
1087 }
1088 break;
1089 }
1090 }
1091 }
1092}
1093
1095 Array<OneD, Array<OneD, NekDouble>> &veldotMF)
1096{
1097 int nq = m_fields[0]->GetNpoints();
1098
1100
1101 Array<OneD, NekDouble> magMF(nq);
1102 for (int j = 0; j < m_shapedim; ++j)
1103 {
1104 veldotMF[j] = Array<OneD, NekDouble>(nq, 0.0);
1105
1106 for (int k = 0; k < m_spacedim; ++k)
1107 {
1108 Vmath::Vvtvp(nq, &m_movingframes[j][k * nq], 1, &m_velocity[k][0],
1109 1, &veldotMF[j][0], 1, &veldotMF[j][0], 1);
1110 }
1111 }
1112}
1113
1114void MMFAdvection::v_EvaluateExactSolution([[maybe_unused]] unsigned int field,
1115 Array<OneD, NekDouble> &outfield,
1116 const NekDouble time)
1117{
1118 switch (m_TestType)
1119 {
1120 case eAdvectionBell:
1121 {
1122 AdvectionBellSphere(outfield);
1123 }
1124 break;
1125
1126 case eTestPlane:
1127 {
1128 Test2Dproblem(time, outfield);
1129 }
1130 break;
1131
1133 {
1134 AdvectionBellPlane(outfield);
1135 }
1136 break;
1137
1138 case eTestCube:
1139 {
1140 Test3Dproblem(time, outfield);
1141 }
1142 break;
1143
1144 default:
1145 break;
1146 }
1147}
1148
1150{
1153 SolverUtils::AddSummaryItem(s, "Rotation Angle", m_RotAngle);
1154}
1155
1156} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
@ 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.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
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:65
void Test2Dproblem(const NekDouble time, 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.
Array< OneD, Array< OneD, NekDouble > > m_veldotMF
Definition: MMFAdvection.h:91
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
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)
void WeakDGDirectionalAdvection(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField)
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
void ComputeNablaCdotVelocity(Array< OneD, NekDouble > &vellc)
void ComputeveldotMF(Array< OneD, Array< OneD, NekDouble > > &veldotMF)
void v_SetInitialConditions(const NekDouble initialtime, bool dumpInitialConditions, const int domain) override
void EvaluateAdvectionVelocity(Array< OneD, Array< OneD, NekDouble > > &velocity)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFAdvection.h:69
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Definition: MMFAdvection.h:83
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Array< OneD, Array< OneD, NekDouble > > m_velocity
Definition: MMFAdvection.h:89
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
void v_DoSolve() override
Solves an unsteady problem.
static std::string className
Name of class.
Definition: MMFAdvection.h:80
Array< OneD, NekDouble > m_vellc
Definition: MMFAdvection.h:92
void AdvectionBellPlane(Array< OneD, NekDouble > &outfield)
Array< OneD, NekDouble > m_traceVn
Definition: MMFAdvection.h:90
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
int m_spacedim
Spatial dimension (>= expansion dim).
NekDouble m_timestep
Time step size.
int m_infosteps
Number of time steps between outputting status information.
NekDouble m_time
Current time of simulation.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
NekDouble m_fintime
Finish time of the simulation.
SOLVER_UTILS_EXPORT int GetNcoeffs()
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.
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:144
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Virtual function for generating summary information.
Definition: MMFSystem.cpp:2463
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:696
SOLVER_UTILS_EXPORT NekDouble AvgInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2292
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:183
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:2399
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:187
SOLVER_UTILS_EXPORT void ComputencdotMF()
Definition: MMFSystem.cpp:317
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:792
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2338
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:186
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:51
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.
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekZeroTol
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
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:47
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::vector< double > z(NPUPPER)
@ 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.hpp:72
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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.hpp:366
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.hpp:180
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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.hpp:220
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285