Nektar++
MMFSWE.cpp
Go to the documentation of this file.
1/////////////////////////////////////////////////////////////////////////////
2//
3// File: MMFSWE.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/predicate.hpp>
39
44
45namespace Nektar
46{
47
48std::string MMFSWE::className =
50 "MMFSWE", MMFSWE::create, "MMFSWE equation.");
51
54 : UnsteadySystem(pSession, pGraph), MMFSystem(pSession, pGraph)
55{
56 m_planeNumber = 0;
57}
58
59/**
60 * @brief Initialisation object for the unsteady linear advection equation.
61 */
62void MMFSWE::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 // Load acceleration of gravity
78 m_session->LoadParameter("Gravity", m_g, 9.81);
79
80 // Add Coriolois effects
81 m_session->LoadParameter("AddCoriolis", m_AddCoriolis, 1);
82
83 // Add Rotation of the sphere along the pole
84 m_session->LoadParameter("AddRotation", m_AddRotation, 1);
85
86 m_session->LoadParameter("AddRossbyDisturbance", m_RossbyDisturbance, 0);
87 m_session->LoadParameter("PurturbedJet", m_PurturbedJet, 0);
88
89 // Define TestType
90 ASSERTL0(m_session->DefinesSolverInfo("TESTTYPE"),
91 "No TESTTYPE defined in session.");
92 std::string TestTypeStr = m_session->GetSolverInfo("TESTTYPE");
93 for (int i = 0; i < (int)SIZE_TestType; ++i)
94 {
95 if (boost::iequals(TestTypeMap[i], TestTypeStr))
96 {
98 break;
99 }
100 }
101
102 // Variable Setting for each test problem
103 NekDouble gms = 9.80616;
104
105 switch (m_TestType)
106 {
107 case eTestSteadyZonal:
108 {
109 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
110 NekDouble rad_earth = 6.37122 * 1000000;
111 NekDouble Omegams;
112
113 // Nondimensionalized coeffs.
114 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
115
116 m_session->LoadParameter("RotationAngle", m_alpha, 0.0);
117 m_session->LoadParameter("u0", m_u0, 2.0 * m_pi / 12.0);
118 m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
119 m_Omega = Omegams * SecondToDay;
120
121 m_session->LoadParameter("H0", m_H0, 2.94 * 10000);
122 m_H0 = m_H0 / (rad_earth * gms);
123
124 m_Hvar = (1.0 / m_g) * (m_Omega * m_u0 + 0.5 * m_u0 * m_u0);
125 }
126 break;
127
129 {
130 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
131 NekDouble rad_earth = 6.37122 * 1000000;
132 NekDouble Omegams;
133
134 // Nondimensionalized coeffs.
135 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
136
137 m_session->LoadParameter("RotationAngle", m_alpha, 0.0);
138 m_session->LoadParameter("u0", m_u0, 2.0 * m_pi / 12.0);
139 m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
140 m_Omega = Omegams * SecondToDay;
141
142 m_H0 = 133681.0 / (rad_earth * gms); // m^2 / s^2
143 m_k2 = 10.0 / (rad_earth * gms); // m^2 / s^2
144 }
145 break;
146
148 {
149 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
150 NekDouble rad_earth = 6.37122 * 1000000;
151 NekDouble Omegams;
152
153 // Nondimensionalized coeffs.
154 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
155
156 m_session->LoadParameter("RotationAngle", m_alpha, 0.0);
157
158 m_session->LoadParameter("u0", m_u0, 20.0);
159 m_u0 = m_u0 * SecondToDay / rad_earth;
160
161 m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
162 m_Omega = Omegams * SecondToDay;
163
164 m_H0 = 5960.0 / rad_earth;
165 m_hs0 = 2000.0 / rad_earth;
166 }
167 break;
168
169 case eTestUnstableJet:
170 {
171 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
172 NekDouble rad_earth = 6.37122 * 1000000;
173 NekDouble Omegams;
174
175 // Nondimensionalized coeffs.
176 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
177
178 m_session->LoadParameter("u0", m_u0, 80.0);
179 m_u0 = m_u0 * SecondToDay / rad_earth;
180
181 m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
182 m_Omega = Omegams * SecondToDay;
183
184 m_session->LoadParameter("H0", m_H0, 10000.0);
185 m_H0 = m_H0 / rad_earth;
186
187 m_uthetamax = 80 * SecondToDay / rad_earth;
188 m_theta0 = m_pi / 7.0;
189 m_theta1 = m_pi / 2.0 - m_theta0;
190 m_en = exp(-4.0 / (m_theta1 - m_theta0) / (m_theta1 - m_theta0));
191 m_hbar = 120.0 / rad_earth;
192
193 std::cout << "m_theta0 = " << m_theta0
194 << ", m_theta1 = " << m_theta1 << ", m_en = " << m_en
195 << ", m_hbar = " << m_hbar << std::endl;
196 }
197 break;
198
199 case eTestRossbyWave:
200 {
201 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
202 NekDouble rad_earth = 6.37122 * 1000000;
203 NekDouble Omegams;
204
205 // Nondimensionalized coeffs.
206 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
207
208 m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
209 m_Omega = Omegams * SecondToDay;
210
211 m_session->LoadParameter("H0", m_H0, 8000.0);
212 m_H0 = m_H0 / rad_earth;
213
214 m_angfreq = 7.848 * 0.000001 * SecondToDay;
215 m_K = 7.848 * 0.000001 * SecondToDay;
216 }
217 break;
218
219 default:
220 break;
221 }
222
223 // TestVorticityComputation
225 {
227 }
228
229 // If explicit it computes RHS and PROJECTION for the time integration
231 {
234 }
235 // Otherwise it gives an error (no implicit integration)
236 else
237 {
238 ASSERTL0(false, "Implicit unsteady Advection not set up.");
239 }
240}
241
243{
244 ASSERTL0(m_intScheme != nullptr, "No time integration scheme.");
245
246 int i, nchk = 1;
247 int nvariables = 0;
248 int nfields = m_fields.size();
249 int nq = m_fields[0]->GetNpoints();
250
251 if (m_intVariables.empty())
252 {
253 for (i = 0; i < nfields; ++i)
254 {
255 m_intVariables.push_back(i);
256 }
257 nvariables = nfields;
258 }
259 else
260 {
261 nvariables = m_intVariables.size();
262 }
263
264 // Set up wrapper to fields data storage.
265 Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
266 Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
267
268 // Order storage to list time-integrated fields first.
269 for (i = 0; i < nvariables; ++i)
270 {
271 fields[i] = m_fields[m_intVariables[i]]->GetPhys();
272 m_fields[m_intVariables[i]]->SetPhysState(false);
273 }
274
275 // Initialise time integration scheme
276 m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
277
278 // Check uniqueness of checkpoint output
279 ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
280 (m_checktime > 0.0 && m_checksteps == 0) ||
281 (m_checktime == 0.0 && m_checksteps > 0),
282 "Only one of IO_CheckTime and IO_CheckSteps "
283 "should be set!");
284
286 bool doCheckTime = false;
287 int step = 0;
288 NekDouble intTime = 0.0;
289 NekDouble cpuTime = 0.0;
290 NekDouble elapsed = 0.0;
291
292 NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
293 Array<OneD, NekDouble> zeta(nq);
294 Array<OneD, Array<OneD, NekDouble>> fieldsprimitive(nvariables);
295 for (int i = 0; i < nvariables; ++i)
296 {
297 fieldsprimitive[i] = Array<OneD, NekDouble>(nq);
298 }
299
301 {
302 timer.Start();
303 fields = m_intScheme->TimeIntegrate(step, m_timestep);
304 timer.Stop();
305
307 elapsed = timer.TimePerTest(1);
308 intTime += elapsed;
309 cpuTime += elapsed;
310
311 // Write out status information
312 if (m_infosteps && !((step + 1) % m_infosteps) &&
313 m_session->GetComm()->GetRank() == 0)
314 {
315 std::cout << "Steps: " << std::setw(8) << std::left << step + 1
316 << " "
317 << "Time: " << std::setw(12) << std::left << m_time;
318
319 std::stringstream ss;
320 ss << cpuTime << "s";
321 std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
322 << std::endl;
323
324 // Printout Mass, Energy, Enstrophy
325 ConservativeToPrimitive(fields, fieldsprimitive);
326
327 // Vorticity zeta
328 ComputeVorticity(fieldsprimitive[1], fieldsprimitive[2], zeta);
329 Vorticity = std::abs(m_fields[0]->Integral(zeta) - m_Vorticity0);
330
331 // Masss = h^*
332 Mass = (ComputeMass(fieldsprimitive[0]) - m_Mass0) / m_Mass0;
333
334 // Energy = 0.5*( h^*(u^2 + v^2) + g ( h^2 - h_s^s ) )
335 Energy = (ComputeEnergy(fieldsprimitive[0], fieldsprimitive[1],
336 fieldsprimitive[2]) -
337 m_Energy0) /
338 m_Energy0;
339
340 // Enstrophy = 0.5/h^* ( \mathbf{k} \cdot (\nabla \times \mathbf{v}
341 // ) + f )^2
342 Enstrophy =
343 (ComputeEnstrophy(fieldsprimitive[0], fieldsprimitive[1],
344 fieldsprimitive[2]) -
345 m_Enstrophy0) /
347
348 std::cout << "dMass = " << std::setw(8) << std::left << Mass << " "
349 << ", dEnergy = " << std::setw(8) << std::left << Energy
350 << " "
351 << ", dEnstrophy = " << std::setw(8) << std::left
352 << Enstrophy << " "
353 << ", dVorticity = " << std::setw(8) << std::left
354 << Vorticity << std::endl
355 << std::endl;
356
357 cpuTime = 0.0;
358 }
359
360 // (h, hu, hv) -> (\eta, u, v)
362
363 // Transform data into coefficient space
364 for (i = 0; i < nvariables; ++i)
365 {
366 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
367 m_fields[m_intVariables[i]]->FwdTransLocalElmt(
368 fields[i], m_fields[m_intVariables[i]]->UpdateCoeffs());
369 m_fields[m_intVariables[i]]->SetPhysState(false);
370 }
371
372 // Write out checkpoint files
373 if ((m_checksteps && step && !((step + 1) % m_checksteps)) ||
374 doCheckTime)
375 {
376 Checkpoint_Output(nchk++);
377 doCheckTime = false;
378
379 // (\eta, u, v) -> (h, hu, hv)
381 }
382
383 // Step advance
384 ++step;
385 }
386
387 // Print out summary statistics
388 if (m_session->GetComm()->GetRank() == 0)
389 {
390 if (m_cflSafetyFactor > 0.0)
391 {
392 std::cout << "CFL safety factor : " << m_cflSafetyFactor
393 << std::endl
394 << "CFL time-step : " << m_timestep << std::endl;
395 }
396
397 if (m_session->GetSolverInfo("Driver") != "SteadyState")
398 {
399 std::cout << "Time-integration : " << intTime << "s" << std::endl;
400 }
401 }
402
403 for (i = 0; i < nvariables; ++i)
404 {
405 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
406 m_fields[m_intVariables[i]]->SetPhysState(true);
407 }
408
409 // (h, hu, hv) -> (\eta, u, v)
411
412 for (i = 0; i < nvariables; ++i)
413 {
414 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
415 m_fields[i]->UpdateCoeffs());
416 }
417}
418
421 [[maybe_unused]] const NekDouble time)
422{
423 int i;
424 int nvariables = inarray.size();
425 int ncoeffs = GetNcoeffs();
426 int nq = GetTotPoints();
427
428 // inarray in physical space
429 Array<OneD, Array<OneD, NekDouble>> physarray(nvariables);
430 Array<OneD, Array<OneD, NekDouble>> modarray(nvariables);
431
432 for (i = 0; i < nvariables; ++i)
433 {
434 physarray[i] = Array<OneD, NekDouble>(nq);
435 modarray[i] = Array<OneD, NekDouble>(ncoeffs);
436 }
437
438 // (h, hu, hv) -> (\eta, u, v)
439 ConservativeToPrimitive(inarray, physarray);
440
441 // Weak Directional Derivative
442 WeakDGSWEDirDeriv(physarray, modarray);
443 AddDivForGradient(physarray, modarray);
444
445 for (i = 0; i < nvariables; ++i)
446 {
447 Vmath::Neg(ncoeffs, modarray[i], 1);
448 }
449
450 // coriolis forcing
451 if (m_AddCoriolis)
452 {
453 AddCoriolis(physarray, modarray);
454 }
455
456 // Bottom Elevation Effect
457 // Add g H \nabla m_depth
458 AddElevationEffect(physarray, modarray);
459
460 // Add terms concerning the rotation of the moving frame
461 if (m_AddRotation)
462 {
463 AddRotation(physarray, modarray);
464 }
465
466 for (i = 0; i < nvariables; ++i)
467 {
468 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
469 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
470 }
471}
472
474 const Array<OneD, Array<OneD, NekDouble>> &InField,
476{
477 int i;
478 int nq = GetNpoints();
479 int ncoeffs = GetNcoeffs();
480 int nTracePointsTot = GetTraceNpoints();
481 int nvariables = m_fields.size();
482
484 Array<OneD, Array<OneD, NekDouble>> physfield(nvariables);
485
486 for (i = 0; i < m_shapedim; ++i)
487 {
488 fluxvector[i] = Array<OneD, NekDouble>(nq);
489 }
490
491 // InField is Primitive
492 for (i = 0; i < nvariables; ++i)
493 {
494 physfield[i] = InField[i];
495 }
496
497 // Get the ith component of the flux vector in (physical space)
498 // fluxvector[0] = component for e^1 cdot \nabla \varphi
499 // fluxvector[1] = component for e^2 cdot \nabla \varphi
500 Array<OneD, NekDouble> fluxtmp(nq);
501 Array<OneD, NekDouble> tmp(ncoeffs);
502
503 // Compute Divergence Components
504 for (i = 0; i < nvariables; ++i)
505 {
506 GetSWEFluxVector(i, physfield, fluxvector);
507
508 OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
509 for (int j = 0; j < m_shapedim; ++j)
510 {
511 // Directional derivation with respect to the j'th moving frame
512 // tmp_j = ( \nabla \phi, fluxvector[j] \mathbf{e}^j )
513 m_fields[i]->IProductWRTDirectionalDerivBase(m_movingframes[j],
514 fluxvector[j], tmp);
515 Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
516 &OutField[i][0], 1);
517 }
518 }
519
520 // Numerical Flux
521 Array<OneD, Array<OneD, NekDouble>> numfluxFwd(nvariables);
522 Array<OneD, Array<OneD, NekDouble>> numfluxBwd(nvariables);
523
524 for (i = 0; i < nvariables; ++i)
525 {
526 numfluxFwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
527 numfluxBwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
528 }
529
530 NumericalSWEFlux(physfield, numfluxFwd, numfluxBwd);
531
532 // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
533 for (i = 0; i < nvariables; ++i)
534 {
535 Vmath::Neg(ncoeffs, OutField[i], 1);
536 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
537 OutField[i]);
538 m_fields[i]->SetPhysState(false);
539 }
540}
541
542// Substract 0.5 * g * H * H / || e^m ||^2 \nalba \cdot e^m
545{
546 // routine works for both primitive and conservative formulations
547 int ncoeffs = outarray[0].size();
548 int nq = physarray[0].size();
549
553 for (int i = 0; i < m_shapedim; ++i)
554 {
555 fluxvector[i] = Array<OneD, NekDouble>(nq);
556 }
557
558 // Get 0.5 g H*H / || e^m ||^2
559 GetSWEFluxVector(3, physarray, fluxvector);
560
561 Array<OneD, NekDouble> tmpc(ncoeffs);
562 for (int j = 0; j < m_shapedim; ++j)
563 {
564 Vmath::Vmul(nq, &fluxvector[0][0], 1, &m_DivMF[j][0], 1, &tmp[0], 1);
565
566 Vmath::Neg(nq, &tmp[0], 1);
567 m_fields[0]->IProductWRTBase(tmp, tmpc);
568
569 Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
570 }
571}
572
574 const int i, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
576{
577 int nq = m_fields[0]->GetTotPoints();
578
579 switch (i)
580 {
581 // flux function for the h equation = [(\eta + d ) u, (\eta + d) v ]
582 case 0:
583 {
584 // h in flux 1
585 Vmath::Vadd(nq, physfield[0], 1, m_depth, 1, flux[1], 1);
586
587 // hu in flux 0
588 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
589
590 // hv in flux 1
591 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
592 }
593 break;
594
595 // flux function for the hu equation = [Hu^2 + 0.5 g H^2, Huv]
596 case 1:
597 {
599
600 // h in tmp
601 Vmath::Vadd(nq, physfield[0], 1, m_depth, 1, tmp, 1);
602
603 // hu in flux 1
604 Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
605
606 // huu in flux 0
607 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
608
609 // hh in tmp
610 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
611
612 // huu + 0.5 g hh in flux 0
613 // Daxpy overwrites flux[0] on exit
614 Blas::Daxpy(nq, 0.5 * m_g, tmp, 1, flux[0], 1);
615
616 // huv in flux 1
617 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
618 }
619 break;
620
621 // flux function for the hv equation = [Huv, Hv^2]
622 case 2:
623 {
625
626 // h in tmp
627 Vmath::Vadd(nq, physfield[0], 1, m_depth, 1, tmp, 1);
628
629 // hv in flux 0
630 Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
631
632 // hvv in flux 1
633 Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
634
635 // huv in flux 0
636 Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
637
638 // hh in tmp
639 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
640
641 // hvv + 0.5 g hh in flux 1
642 Blas::Daxpy(nq, 0.5 * m_g, tmp, 1, flux[1], 1);
643 }
644 break;
645
646 // flux function 0.5 g h * h
647 case 3:
648 {
650 flux[0] = Array<OneD, NekDouble>(nq, 0.0);
651
652 // h in tmp
653 Vmath::Vadd(nq, physfield[0], 1, m_depth, 1, h, 1);
654
655 // hh in tmp
656 Vmath::Vmul(nq, h, 1, h, 1, h, 1);
657
658 // 0.5 g hh in flux 0
659 Blas::Daxpy(nq, 0.5 * m_g, h, 1, flux[0], 1);
660 }
661 break;
662
663 default:
664 break;
665 }
666}
667
669 Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
670 Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
671{
672
673 int i, k;
674 int nTraceNumPoints = GetTraceTotPoints();
675 int nvariables = 3; // only the dependent variables
676
677 // get temporary arrays
678 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
679 Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
680 Array<OneD, NekDouble> DepthFwd(nTraceNumPoints);
681 Array<OneD, NekDouble> DepthBwd(nTraceNumPoints);
682
683 for (i = 0; i < nvariables; ++i)
684 {
685 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
686 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
687 }
688
689 // get the physical values at the trace
690 for (i = 0; i < nvariables; ++i)
691 {
692 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
693
694 // Copy Fwd to Bwd at boundaries
696 }
697 m_fields[0]->GetFwdBwdTracePhys(m_depth, DepthFwd, DepthBwd);
698 CopyBoundaryTrace(DepthFwd, DepthBwd, SolverUtils::eFwdEQBwd);
699
700 // note that we are using the same depth - i.e. the depth is assumed
701 // continuous...
702 switch (m_upwindType)
703 {
705 {
706 // rotate the values to the normal direction
707 NekDouble tmpX, tmpY;
708
709 // Fwd[1] = hu^+ = ( h^1^+ (e^1 \cdot n) + h^2^+ ( e^2 \cdot n) )
710 // Fwd[2] = hv^+ = ( h^1^+ (e^1 \cdot n^{\perp}) + h^2^+ ( e^2 \cdot
711 // n^{\perp}) )
712 // Bwd[1] = hu^- = ( h^1^- (e^1 \cdot n) + h^2^- ( e^2 \cdot n) )
713 // Bwd[2] = hv^- = ( h^1^- (e^1 \cdot n^{\perp}) + h^2^- ( e^2 \cdot
714 // n^{\perp}) )
715
716 for (k = 0; k < nTraceNumPoints; ++k)
717 {
718 tmpX = Fwd[1][k] * m_ncdotMFFwd[0][k] +
719 Fwd[2][k] * m_ncdotMFFwd[1][k];
720 tmpY = Fwd[1][k] * m_nperpcdotMFFwd[0][k] +
721 Fwd[2][k] * m_nperpcdotMFFwd[1][k];
722 Fwd[1][k] = tmpX;
723 Fwd[2][k] = tmpY;
724
725 tmpX = Bwd[1][k] * m_ncdotMFBwd[0][k] +
726 Bwd[2][k] * m_ncdotMFBwd[1][k];
727 tmpY = Bwd[1][k] * m_nperpcdotMFBwd[0][k] +
728 Bwd[2][k] * m_nperpcdotMFBwd[1][k];
729 Bwd[1][k] = tmpX;
730 Bwd[2][k] = tmpY;
731 }
732
733 // Solve the Riemann problem
734 NekDouble denomFwd, denomBwd;
735 Array<OneD, NekDouble> numfluxF(nvariables);
736 Array<OneD, NekDouble> numfluxB(nvariables);
737
738 NekDouble eF1n, eF2n, eB1n, eB2n;
739 NekDouble eF1t, eF2t, eB1t, eB2t;
740 for (k = 0; k < nTraceNumPoints; ++k)
741 {
742 RiemannSolverHLLC(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
743 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
744 Bwd[2][k], numfluxF, numfluxB);
745
746 // uflux = 1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
747 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( e^2 \cdot n^{\perp} huflux
748 // - e^2 \cdot n hvflux
749 // vflux = -1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
750 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( -e^1 \cdot n^{\perp}
751 // huflux + e^1 \cdot n hvflux
752 eF1n = m_ncdotMFFwd[0][k];
753 eF2n = m_ncdotMFFwd[1][k];
754 eB1n = m_ncdotMFBwd[0][k];
755 eB2n = m_ncdotMFBwd[1][k];
756
757 eF1t = m_nperpcdotMFFwd[0][k];
758 eF2t = m_nperpcdotMFFwd[1][k];
759 eB1t = m_nperpcdotMFBwd[0][k];
760 eB2t = m_nperpcdotMFBwd[1][k];
761
762 denomFwd = eF1n * eF2t - eF2n * eF1t;
763 denomBwd = eB1n * eB2t - eB2n * eB1t;
764
765 numfluxFwd[0][k] = numfluxF[0];
766 numfluxFwd[1][k] = (1.0 / denomFwd) *
767 (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
768 numfluxFwd[2][k] =
769 (1.0 / denomFwd) *
770 (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
771
772 numfluxBwd[0][k] = 1.0 * numfluxB[0];
773 numfluxBwd[1][k] = (1.0 / denomBwd) *
774 (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
775 numfluxBwd[2][k] =
776 (1.0 / denomBwd) *
777 (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
778 }
779 }
780 break;
781
785 {
786 Array<OneD, NekDouble> numfluxF(nvariables * (m_shapedim + 1));
787 Array<OneD, NekDouble> numfluxB(nvariables * (m_shapedim + 1));
788
789 for (k = 0; k < nTraceNumPoints; ++k)
790 {
792 {
793 AverageFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
794 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
795 Bwd[2][k], numfluxF, numfluxB);
796 }
797
799 {
800 LaxFriedrichFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
801 Fwd[2][k], Bwd[0][k] + DepthFwd[k],
802 Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
803 }
804
806 {
807 RusanovFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
808 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
809 Bwd[2][k], numfluxF, numfluxB);
810 }
811
812 int indx;
813 NekDouble tmpF0, tmpF1, tmpB0, tmpB1;
814 for (i = 0; i < nvariables; ++i)
815 {
816 indx = i * (m_shapedim + 1);
817
818 tmpF0 = numfluxF[indx] * m_ncdotMFFwd[0][k];
819 tmpF1 = numfluxF[indx + 1] * m_ncdotMFFwd[1][k];
820
821 tmpB0 = numfluxB[indx] * m_ncdotMFBwd[0][k];
822 tmpB1 = numfluxB[indx + 1] * m_ncdotMFBwd[1][k];
823
824 numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
825 numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
826 }
827 }
828 }
829 break;
830
831 default:
832 {
833 ASSERTL0(false, "populate switch statement for upwind flux");
834 }
835 break;
836 }
837}
838
839void MMFSWE::RiemannSolverHLLC([[maybe_unused]] const int index, NekDouble hL,
840 NekDouble uL, NekDouble vL, NekDouble hR,
841 NekDouble uR, NekDouble vR,
842 Array<OneD, NekDouble> &numfluxF,
843 Array<OneD, NekDouble> &numfluxB)
844{
845 NekDouble g = m_g;
846
847 NekDouble cL = sqrt(g * hL);
848 NekDouble cR = sqrt(g * hR);
849
850 NekDouble uRF, vRF, uLB, vLB;
851 NekDouble hstarF, hstarB;
852
853 // Temporary assignments
854 uRF = uR;
855 vRF = vR;
856 uLB = uL;
857 vLB = vL;
858
859 // the two-rarefaction wave assumption
860 hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
861 hstarF *= hstarF;
862 hstarF *= (1.0 / g);
863
864 hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
865 hstarB *= hstarB;
866 hstarB *= (1.0 / g);
867
868 NekDouble hfluxF, hufluxF, hvfluxF;
869 NekDouble hfluxB, hufluxB, hvfluxB;
870 Computehhuhvflux(hL, uL, vL, hR, uRF, vRF, hstarF, hfluxF, hufluxF,
871 hvfluxF);
872 Computehhuhvflux(hL, uLB, vLB, hR, uR, vR, hstarB, hfluxB, hufluxB,
873 hvfluxB);
874
875 numfluxF[0] = hfluxF;
876 numfluxF[1] = hufluxF;
877 numfluxF[2] = hvfluxF;
878
879 numfluxB[0] = hfluxB;
880 numfluxB[1] = hufluxB;
881 numfluxB[2] = hvfluxB;
882}
883
885 NekDouble hR, NekDouble uR, NekDouble vR,
886 NekDouble hstar, NekDouble &hflux,
887 NekDouble &huflux, NekDouble &hvflux)
888{
889 NekDouble g = m_g;
890
891 NekDouble hC, huC, hvC, SL, SR, Sstar;
892 NekDouble cL = sqrt(g * hL);
893 NekDouble cR = sqrt(g * hR);
894
895 // Compute SL
896 if (hstar > hL)
897 {
898 SL = uL - cL * sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
899 }
900 else
901 {
902 SL = uL - cL;
903 }
904
905 // Compute SR
906 if (hstar > hR)
907 {
908 SR = uR + cR * sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
909 }
910 else
911 {
912 SR = uR + cR;
913 }
914
915 if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
916 {
917 Sstar = 0.0;
918 }
919 else
920 {
921 Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
922 (hR * (uR - SR) - hL * (uL - SL));
923 }
924
925 if (SL >= 0)
926 {
927 hflux = hL * uL;
928 huflux = uL * uL * hL + 0.5 * g * hL * hL;
929 hvflux = hL * uL * vL;
930 }
931 else if (SR <= 0)
932 {
933 hflux = hR * uR;
934 huflux = uR * uR * hR + 0.5 * g * hR * hR;
935 hvflux = hR * uR * vR;
936 }
937 else
938 {
939 if ((SL < 0) && (Sstar >= 0))
940 {
941 hC = hL * ((SL - uL) / (SL - Sstar));
942 huC = hC * Sstar;
943 hvC = hC * vL;
944
945 hflux = hL * uL + SL * (hC - hL);
946 huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
947 hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
948 }
949 else
950 {
951 hC = hR * ((SR - uR) / (SR - Sstar));
952 huC = hC * Sstar;
953 hvC = hC * vR;
954
955 hflux = hR * uR + SR * (hC - hR);
956 huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
957 hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
958 }
959 }
960}
961
962void MMFSWE::AverageFlux(const int index, NekDouble hL, NekDouble uL,
963 NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR,
964 Array<OneD, NekDouble> &numfluxF,
965 Array<OneD, NekDouble> &numfluxB)
966{
967 NekDouble MageF1, MageF2, MageB1, MageB2;
968 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
969 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
970
971 NekDouble g = m_g;
972 NekDouble uRF, vRF, uLB, vLB;
973
974 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
975 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
976
977 // uRF = uR component in moving frames e^{Fwd}
978 // vRF = vR component in moving frames e^{Fwd}
979 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
980 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
981
982 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
983 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
984 numfluxF[2] = 0.0;
985
986 numfluxF[3] =
987 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
988 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
989 numfluxF[5] = 0.0;
990
991 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
992 numfluxF[7] =
993 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
994 numfluxF[8] = 0.0;
995
996 // uLB = uL component in moving frames e^{Bwd}
997 // vLB = vL component in moving frames e^{Bwd}
998 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
999 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1000
1001 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1002 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1003 numfluxB[2] = 0.0;
1004
1005 numfluxB[3] =
1006 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1007 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1008 numfluxB[5] = 0.0;
1009
1010 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1011 numfluxB[7] =
1012 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1013 numfluxB[8] = 0.0;
1014}
1015
1016void MMFSWE::LaxFriedrichFlux(const int index, NekDouble hL, NekDouble uL,
1017 NekDouble vL, NekDouble hR, NekDouble uR,
1018 NekDouble vR, Array<OneD, NekDouble> &numfluxF,
1019 Array<OneD, NekDouble> &numfluxB)
1020{
1021 int nvariables = 3;
1022 NekDouble MageF1, MageF2, MageB1, MageB2;
1023 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1024 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1025
1026 NekDouble g = m_g;
1027 NekDouble uRF, vRF, uLB, vLB;
1028 NekDouble velL, velR, lambdaF, lambdaB;
1029
1030 Array<OneD, NekDouble> EigF(nvariables);
1031 Array<OneD, NekDouble> EigB(nvariables);
1032
1033 // Compute Magnitude and Dot product of moving frames for the index
1034 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1035 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1036
1037 // Get the velocity in the normal to the edge
1038 velL = uL * m_ncdotMFFwd[0][index] + vL * m_ncdotMFFwd[1][index];
1039 velR = -1.0 * (uR * m_ncdotMFBwd[0][index] + vR * m_ncdotMFBwd[1][index]);
1040
1041 EigF[0] = velL - sqrt(g * hL);
1042 EigF[1] = velL;
1043 EigF[2] = velL + sqrt(g * hL);
1044
1045 EigB[0] = velR - sqrt(g * hR);
1046 EigB[1] = velR;
1047 EigB[2] = velR + sqrt(g * hR);
1048
1049 lambdaF = Vmath::Vamax(nvariables, EigF, 1);
1050 lambdaB = Vmath::Vamax(nvariables, EigB, 1);
1051
1052 // uRF = uR component in moving frames e^{Fwd}
1053 // vRF = vR component in moving frames e^{Fwd}
1054 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1055 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1056
1057 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1058 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1059 numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1060
1061 numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1062 0.5 * g * (hL * hL + hR * hR));
1063 numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1064 numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1065
1066 numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1067 numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1068 0.5 * g * (hL * hL + hR * hR));
1069 numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1070
1071 // uLB = uL component in moving frames e^{Bwd}
1072 // vLB = vL component in moving frames e^{Bwd}
1073 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1074 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1075
1076 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1077 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1078 numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1079
1080 numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1081 0.5 * g * (hR * hR + hL * hL));
1082 numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1083 numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1084
1085 numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1086 numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1087 0.5 * g * (hR * hR + hL * hL));
1088 numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1089}
1090
1091void MMFSWE::RusanovFlux(const int index, NekDouble hL, NekDouble uL,
1092 NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR,
1093 Array<OneD, NekDouble> &numfluxF,
1094 Array<OneD, NekDouble> &numfluxB)
1095{
1096 int nvariables = 3;
1097 NekDouble MageF1, MageF2, MageB1, MageB2;
1098 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1099 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1100
1101 NekDouble g = m_g;
1102 NekDouble uRF, vRF, uLB, vLB;
1103 NekDouble velL, velR;
1104
1105 Array<OneD, NekDouble> EigF(nvariables);
1106 Array<OneD, NekDouble> EigB(nvariables);
1107
1108 // Compute Magnitude and Dot product of moving frames for the index
1109 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1110 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1111
1112 // Get the velocity in the normal to the edge
1113 velL = uL * m_ncdotMFFwd[0][index] + vL * m_ncdotMFFwd[1][index];
1114 velR = -1.0 * (uR * m_ncdotMFBwd[0][index] + vR * m_ncdotMFBwd[1][index]);
1115
1116 NekDouble SL, SR;
1117 SL = fabs(velL) + sqrt(g * hL);
1118 SR = fabs(velR) + sqrt(g * hR);
1119
1120 NekDouble S;
1121 if (SL > SR)
1122 {
1123 S = SL;
1124 }
1125 else
1126 {
1127 S = SR;
1128 }
1129
1130 // uRF = uR component in moving frames e^{Fwd}
1131 // vRF = vR component in moving frames e^{Fwd}
1132 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1133 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1134
1135 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1136 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1137 numfluxF[2] = 0.5 * S * (hL - hR);
1138
1139 numfluxF[3] =
1140 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1141 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1142 numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1143
1144 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1145 numfluxF[7] =
1146 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1147 numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1148
1149 // uLB = uL component in moving frames e^{Bwd}
1150 // vLB = vL component in moving frames e^{Bwd}
1151 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1152 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1153
1154 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1155 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1156 numfluxB[2] = 0.5 * S * (hL - hR);
1157
1158 numfluxB[3] =
1159 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1160 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1161 numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1162
1163 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1164 numfluxB[7] =
1165 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1166 numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1167}
1168
1169void MMFSWE::ComputeMagAndDot(const int index, NekDouble &MageF1,
1170 NekDouble &MageF2, NekDouble &MageB1,
1171 NekDouble &MageB2, NekDouble &eF1_cdot_eB1,
1172 NekDouble &eF1_cdot_eB2, NekDouble &eF2_cdot_eB1,
1173 NekDouble &eF2_cdot_eB2)
1174{
1175 NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1176 NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1177
1178 MF1x = m_MFtraceFwd[0][0][index];
1179 MF1y = m_MFtraceFwd[0][1][index];
1180 MF1z = m_MFtraceFwd[0][2][index];
1181
1182 MF2x = m_MFtraceFwd[1][0][index];
1183 MF2y = m_MFtraceFwd[1][1][index];
1184 MF2z = m_MFtraceFwd[1][2][index];
1185
1186 MB1x = m_MFtraceBwd[0][0][index];
1187 MB1y = m_MFtraceBwd[0][1][index];
1188 MB1z = m_MFtraceBwd[0][2][index];
1189
1190 MB2x = m_MFtraceBwd[1][0][index];
1191 MB2y = m_MFtraceBwd[1][1][index];
1192 MB2z = m_MFtraceBwd[1][2][index];
1193
1194 // MFtrace = MFtrace [ j*spacedim + k ], j = shape, k = sapce
1195 MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1196 MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1197 MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1198 MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1199
1200 eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1201 eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1202 eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1203 eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1204}
1205
1206// Add Coriolis factors
1208 Array<OneD, Array<OneD, NekDouble>> &outarray)
1209{
1210 int ncoeffs = outarray[0].size();
1211 int nq = physarray[0].size();
1212
1214 Array<OneD, NekDouble> tmp(nq);
1215 Array<OneD, NekDouble> tmpc(ncoeffs);
1216
1217 // physarray is primitive
1218 // conservative formulation compute h
1219 // h = \eta + d
1220 Vmath::Vadd(nq, physarray[0], 1, m_depth, 1, h, 1);
1221
1222 int indx = 0;
1223 for (int j = 0; j < m_shapedim; ++j)
1224 {
1225 if (j == 0)
1226 {
1227 indx = 2;
1228 }
1229
1230 else if (j == 1)
1231 {
1232 indx = 1;
1233 }
1234
1235 // add to hu equation
1236 Vmath::Vmul(nq, m_coriolis, 1, physarray[indx], 1, tmp, 1);
1237 Vmath::Vmul(nq, h, 1, tmp, 1, tmp, 1);
1238
1239 if (j == 1)
1240 {
1241 Vmath::Neg(nq, tmp, 1);
1242 }
1243
1244 // N \cdot (e^1 \times e^2 )
1245 // Vmath::Vmul(nq, &m_MF1crossMF2dotSN[0], 1, &tmp[0], 1, &tmp[0], 1);
1246 m_fields[0]->IProductWRTBase(tmp, tmpc);
1247 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1248 }
1249}
1250
1251// Compuate g H \nabla m_depth
1253 Array<OneD, Array<OneD, NekDouble>> &outarray)
1254{
1255 int ncoeffs = outarray[0].size();
1256 int nq = physarray[0].size();
1257
1259 Array<OneD, NekDouble> tmp(nq);
1260 Array<OneD, NekDouble> tmpc(ncoeffs);
1261
1262 // physarray is primitive
1263 // conservative formulation compute h
1264 // h = \eta + d
1265 Vmath::Vadd(nq, physarray[0], 1, m_depth, 1, h, 1);
1266
1267 for (int j = 0; j < m_shapedim; ++j)
1268 {
1269 Vmath::Vmul(nq, &h[0], 1, &m_Derivdepth[j][0], 1, &tmp[0], 1);
1270 Vmath::Smul(nq, m_g, tmp, 1, tmp, 1);
1271
1272 m_fields[0]->IProductWRTBase(tmp, tmpc);
1273
1274 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1275 }
1276}
1277
1278// =================================================
1279// Add rotational factors
1280// =================================================
1282 Array<OneD, Array<OneD, NekDouble>> &outarray)
1283{
1284 // routine works for both primitive and conservative formulations
1285 int ncoeffs = outarray[0].size();
1286 int nq = physarray[0].size();
1287
1288 // Compute h
1290 Vmath::Vadd(nq, &physarray[0][0], 1, &m_depth[0], 1, &h[0], 1);
1291
1292 Array<OneD, NekDouble> de0dt_cdot_e0;
1293 Array<OneD, NekDouble> de0dt_cdot_e1;
1294 Array<OneD, NekDouble> de1dt_cdot_e0;
1295 Array<OneD, NekDouble> de1dt_cdot_e1;
1296 Compute_demdt_cdot_ek(0, 0, physarray, de0dt_cdot_e0);
1297 Compute_demdt_cdot_ek(1, 0, physarray, de1dt_cdot_e0);
1298 Compute_demdt_cdot_ek(0, 1, physarray, de0dt_cdot_e1);
1299 Compute_demdt_cdot_ek(1, 1, physarray, de1dt_cdot_e1);
1300
1301 Array<OneD, NekDouble> Rott1(nq);
1302 Array<OneD, NekDouble> Rott2(nq);
1303 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1304 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1305 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1306 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1307
1308 // Multiply H and \partial \phi / \partial t which is assumed to be u_{\phi}
1309 Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1310 Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1311
1312 Vmath::Neg(nq, Rott1, 1);
1313 Vmath::Neg(nq, Rott2, 1);
1314
1315 Array<OneD, NekDouble> tmpc1(ncoeffs);
1316 Array<OneD, NekDouble> tmpc2(ncoeffs);
1317 m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1318 m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1319
1320 Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1321 Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1322}
1323
1324// =====================================================================
1325// Compute \frac{d e^{indm}}{dt} \cdot e^{indk} / \| e^{indk} \|
1326// Equivalently, [ \frac{d e^{indm}}{d \xi_1} \frac{ d \xi_1}{d \phi} + \frac{d
1327// e^{indm}}{d \xi_2} \frac{ d \xi_2}{d \phi} ] \frac{d \phi}{dt}
1329 const int indm, const int indk,
1330 const Array<OneD, const Array<OneD, NekDouble>> &physarray,
1331 Array<OneD, NekDouble> &outarray)
1332{
1333 int j, k;
1334 int nq = m_fields[0]->GetNpoints();
1335
1336 Array<OneD, NekDouble> tmp(nq, 0.0);
1337 Array<OneD, NekDouble> tmpderiv(nq);
1338
1339 outarray = Array<OneD, NekDouble>(nq, 0.0);
1340 for (j = 0; j < m_shapedim; ++j)
1341 {
1342 for (k = 0; k < m_spacedim; ++k)
1343 {
1344 // Compute d e^m / d \xi_1 and d e^m / d \xi_2
1345 Vmath::Vcopy(nq, &m_movingframes[indm][k * nq], 1, &tmp[0], 1);
1346 m_fields[0]->PhysDirectionalDeriv(m_movingframes[j], tmp, tmpderiv);
1347
1348 Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1349 &tmpderiv[0], 1);
1350
1351 Vmath::Vvtvp(nq, &tmpderiv[0], 1, &m_movingframes[indk][k * nq], 1,
1352 &outarray[0], 1, &outarray[0], 1);
1353 }
1354 }
1355}
1356
1357/**
1358 * @brief Compute the projection for the linear advection equation.
1359 *
1360 * @param inarray Given fields.
1361 * @param outarray Calculated solution.
1362 * @param time Time.
1363 */
1365 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1366 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
1367{
1368 switch (m_projectionType)
1369 {
1371 {
1372 ConservativeToPrimitive(inarray, outarray);
1373 SetBoundaryConditions(outarray, time);
1374 PrimitiveToConservative(outarray, outarray);
1375 }
1376 break;
1377 default:
1378 ASSERTL0(false, "Unknown projection scheme");
1379 break;
1380 }
1381}
1382
1384 NekDouble time)
1385{
1386
1387 int nvariables = m_fields.size();
1388 int cnt = 0;
1389
1390 // loop over Boundary Regions
1391 for (int n = 0; n < m_fields[0]->GetBndConditions().size(); ++n)
1392 {
1393
1394 // Zonal Boundary Condition
1395 if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() == "eMG")
1396 {
1397 if (m_expdim == 1)
1398 {
1399 ASSERTL0(false, "Illegal dimension");
1400 }
1401 else if (m_expdim == 2)
1402 {
1403 // ZonalBoundary2D(n,cnt,inarray);
1404 }
1405 }
1406
1407 // Wall Boundary Condition
1408 if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() == "eWall")
1409 {
1410 if (m_expdim == 1)
1411 {
1412 ASSERTL0(false, "Illegal dimension");
1413 }
1414 else if (m_expdim == 2)
1415 {
1416 WallBoundary2D(n, cnt, inarray);
1417 }
1418 }
1419
1420 // Time Dependent Boundary Condition (specified in meshfile)
1421 if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
1422 "eTimeDependent")
1423 {
1424 for (int i = 0; i < nvariables; ++i)
1425 {
1426 m_fields[i]->EvaluateBoundaryConditions(time);
1427 }
1428 }
1429 cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1430 }
1431}
1432
1433void MMFSWE::WallBoundary2D(int bcRegion, int cnt,
1434 Array<OneD, Array<OneD, NekDouble>> &physarray)
1435{
1436
1437 int i;
1438 int nTraceNumPoints = GetTraceTotPoints();
1439 int nvariables = physarray.size();
1440
1441 // get physical values of the forward trace
1442 Array<OneD, Array<OneD, NekDouble>> Fwd0(nvariables);
1443 for (i = 0; i < nvariables; ++i)
1444 {
1445 Fwd0[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1446 m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1447 }
1448
1449 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
1450 for (i = 0; i < nvariables; ++i)
1451 {
1452 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1453 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1454 }
1455
1456 // Adjust the physical values of the trace to take
1457 // user defined boundaries into account
1458 int e, id1, id2, npts;
1459
1460 for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
1461 ++e)
1462 {
1463 npts = m_fields[0]
1464 ->GetBndCondExpansions()[bcRegion]
1465 ->GetExp(e)
1466 ->GetTotPoints();
1467 id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1468 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
1469 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1470
1471 switch (m_expdim)
1472 {
1473 case 1:
1474 {
1475 // negate the forward flux
1476 Vmath::Neg(npts, &Fwd[1][id2], 1);
1477 }
1478 break;
1479 case 2:
1480 {
1481 Array<OneD, NekDouble> tmp_n(npts);
1482 Array<OneD, NekDouble> tmp_t(npts);
1483
1484 Vmath::Vmul(npts, &Fwd[1][id2], 1, &m_ncdotMFFwd[0][id2], 1,
1485 &tmp_n[0], 1);
1486 Vmath::Vvtvp(npts, &Fwd[2][id2], 1, &m_ncdotMFFwd[1][id2], 1,
1487 &tmp_n[0], 1, &tmp_n[0], 1);
1488
1489 Vmath::Vmul(npts, &Fwd[1][id2], 1, &m_nperpcdotMFFwd[0][id2], 1,
1490 &tmp_t[0], 1);
1491 Vmath::Vvtvp(npts, &Fwd[2][id2], 1, &m_nperpcdotMFFwd[1][id2],
1492 1, &tmp_t[0], 1, &tmp_t[0], 1);
1493
1494 // negate the normal flux
1495 Vmath::Neg(npts, tmp_n, 1);
1496
1497 Array<OneD, NekDouble> denom(npts);
1498 Array<OneD, NekDouble> tmp_u(npts);
1499 Array<OneD, NekDouble> tmp_v(npts);
1500
1501 // denom = (e^1 \cdot n ) (e^2 \cdot t) - (e^2 \cdot n ) (e^1
1502 // \cdot t)
1503 Vmath::Vmul(npts, &m_ncdotMFFwd[1][id2], 1,
1504 &m_nperpcdotMFFwd[0][id2], 1, &denom[0], 1);
1505 Vmath::Vvtvm(npts, &m_ncdotMFFwd[0][id2], 1,
1506 &m_nperpcdotMFFwd[1][id2], 1, &denom[0], 1,
1507 &denom[0], 1);
1508
1509 Vmath::Vmul(npts, &m_ncdotMFFwd[1][id2], 1, &tmp_t[0], 1,
1510 &tmp_u[0], 1);
1511 Vmath::Vvtvm(npts, &m_nperpcdotMFFwd[1][id2], 1, &tmp_n[0], 1,
1512 &tmp_u[0], 1, &tmp_u[0], 1);
1513 Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1514
1515 Vmath::Vcopy(npts, &tmp_u[0], 1, &Fwd[1][id2], 1);
1516
1517 Vmath::Vmul(npts, &m_nperpcdotMFFwd[0][id2], 1, &tmp_n[0], 1,
1518 &tmp_v[0], 1);
1519 Vmath::Vvtvm(npts, &m_ncdotMFFwd[0][id2], 1, &tmp_t[0], 1,
1520 &tmp_v[0], 1, &tmp_v[0], 1);
1521 Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1522
1523 Vmath::Vcopy(npts, &tmp_v[0], 1, &Fwd[2][id2], 1);
1524 }
1525 break;
1526
1527 default:
1528 ASSERTL0(false, "Illegal expansion dimension");
1529 }
1530
1531 // copy boundary adjusted values into the boundary expansion
1532 for (i = 0; i < nvariables; ++i)
1533 {
1534 Vmath::Vcopy(npts, &Fwd[i][id2], 1,
1535 &(m_fields[i]
1536 ->GetBndCondExpansions()[bcRegion]
1537 ->UpdatePhys())[id1],
1538 1);
1539 }
1540 }
1541}
1542
1543/**
1544 *
1545 */
1546void MMFSWE::v_DoInitialise(bool dumpInitialConditions)
1547{
1548 // Compute m_depth and m_Derivdepth
1551 SetInitialConditions(0.0, dumpInitialConditions);
1553
1554 // transfer the initial conditions to modal values
1555 for (int i = 0; i < m_fields.size(); ++i)
1556 {
1557 m_fields[i]->SetPhysState(true);
1558 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1559 m_fields[i]->UpdateCoeffs());
1560 }
1561}
1562
1564{
1565 int nq = GetTotPoints();
1566
1568
1569 switch (m_TestType)
1570 {
1571 case eTestPlane:
1572 {
1573 Vmath::Fill(nq, 1.0, m_depth, 1);
1574 }
1575 break;
1576
1577 case eTestUnsteadyZonal:
1578 {
1579 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1583
1584 m_fields[0]->GetCoords(x, y, z);
1585
1586 NekDouble x0j, x1j, x2j;
1587 NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1588
1589 for (int j = 0; j < nq; ++j)
1590 {
1591 x0j = x[j];
1592 x1j = y[j];
1593 x2j = z[j];
1594
1595 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1596 sin_theta, cos_theta);
1597 m_depth[j] =
1598 m_H0 - m_k2 -
1599 (0.5 / m_g) * (m_Omega * sin_theta) * (m_Omega * sin_theta);
1600 }
1601 }
1602 break;
1603
1605 {
1606 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1610
1611 m_fields[0]->GetCoords(x, y, z);
1612
1613 int indx = 0;
1614 NekDouble x0j, x1j, x2j, dist2;
1615 NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1616 NekDouble hRad, phic, thetac;
1617 NekDouble Tol = 0.000001;
1618
1619 hRad = m_pi / 9.0;
1620 phic = -m_pi / 2.0;
1621 thetac = m_pi / 6.0;
1622
1623 for (int j = 0; j < nq; ++j)
1624 {
1625 x0j = x[j];
1626 x1j = y[j];
1627 x2j = z[j];
1628
1629 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1630 sin_theta, cos_theta);
1631
1632 if ((std::abs(sin(phic) - sin_varphi) +
1633 std::abs(sin(thetac) - sin_theta)) < Tol)
1634 {
1635 std::cout << "A point " << j
1636 << " is coincient with the singularity "
1637 << std::endl;
1638 indx = 1;
1639 }
1640
1641 phi = atan2(sin_varphi, cos_varphi);
1642 theta = atan2(sin_theta, cos_theta);
1643
1644 // Compute r
1645 dist2 = (phi - phic) * (phi - phic) +
1646 (theta - thetac) * (theta - thetac);
1647
1648 if (dist2 > hRad * hRad)
1649 {
1650 dist2 = hRad * hRad;
1651 }
1652
1653 m_depth[j] = m_H0 - m_hs0 * (1.0 - sqrt(dist2) / hRad);
1654 }
1655
1656 if (!indx)
1657 {
1658 std::cout << "No point is coincident with the singularity point"
1659 << std::endl;
1660 }
1661 }
1662 break;
1663
1664 case eTestUnstableJet:
1665 {
1666 for (int j = 0; j < nq; ++j)
1667 {
1668 m_depth[j] = m_H0;
1669 }
1670 }
1671 break;
1672
1673 case eTestSteadyZonal:
1674 case eTestRossbyWave:
1675 {
1676 Vmath::Zero(nq, m_depth, 1);
1677 }
1678 break;
1679
1680 default:
1681 {
1682 Vmath::Zero(nq, m_depth, 1);
1683 }
1684 break;
1685 }
1686
1687 // Comptue \nabla m_depth \cdot e^i
1689
1690 for (int j = 0; j < m_shapedim; j++)
1691 {
1693 m_fields[0]->PhysDirectionalDeriv(m_movingframes[j], m_depth,
1694 m_Derivdepth[j]);
1695 }
1696
1697 std::cout << "Water Depth (m_depth) was generated with mag = "
1698 << AvgAbsInt(m_depth) << " with max. deriv = ( "
1699 << Vmath::Vamax(nq, m_Derivdepth[0], 1) << " , "
1700 << Vmath::Vamax(nq, m_Derivdepth[1], 1) << " ) " << std::endl;
1701}
1702
1704{
1705 switch (m_TestType)
1706 {
1707 case eTestPlane:
1708 {
1709 GetFunction("Coriolis")->Evaluate("f", m_coriolis);
1710 }
1711 break;
1712
1713 case eTestSteadyZonal:
1714 {
1716 }
1717 break;
1718
1719 case eTestUnsteadyZonal:
1721 case eTestUnstableJet:
1722 case eTestRossbyWave:
1723 {
1725 }
1726 break;
1727
1728 default:
1729 break;
1730 }
1731}
1732
1734{
1735 int nq = GetTotPoints();
1736 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1737
1741
1742 m_fields[0]->GetCoords(x, y, z);
1743
1744 NekDouble x0j, x1j, x2j;
1745 NekDouble tmp;
1746
1747 outarray = Array<OneD, NekDouble>(nq, 0.0);
1748 for (int j = 0; j < nq; ++j)
1749 {
1750 x0j = x[j];
1751 x1j = y[j];
1752 x2j = z[j];
1753
1754 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1755 cos_theta);
1756
1757 // H = 2 \Omega *(- \cos \phi \cos \theta \sin \alpha + \sin \theta \cos
1758 // \alpha )
1759 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
1760 sin_theta * cos(m_alpha);
1761 outarray[j] = 2.0 * m_Omega * tmp;
1762 }
1763}
1764
1766{
1767 int nq = GetTotPoints();
1768
1769 NekDouble x0j, x1j, x2j;
1770 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1771
1775
1776 m_fields[0]->GetCoords(x, y, z);
1777
1778 outarray = Array<OneD, NekDouble>(nq, 0.0);
1779 for (int j = 0; j < nq; ++j)
1780 {
1781 x0j = x[j];
1782 x1j = y[j];
1783 x2j = z[j];
1784
1785 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1786 cos_theta);
1787
1788 outarray[j] = 2.0 * m_Omega * sin_theta;
1789 }
1790}
1791
1793 bool dumpInitialConditions,
1794 [[maybe_unused]] const int domain)
1795{
1796 int nq = GetTotPoints();
1797
1798 switch (m_TestType)
1799 {
1800 case eTestPlane:
1801 {
1802 Array<OneD, NekDouble> eta0(nq);
1805
1806 TestSWE2Dproblem(initialtime, 0, eta0);
1807 m_fields[0]->SetPhys(eta0);
1808
1809 TestSWE2Dproblem(initialtime, 1, u0);
1810 m_fields[1]->SetPhys(u0);
1811
1812 TestSWE2Dproblem(initialtime, 2, v0);
1813 m_fields[2]->SetPhys(v0);
1814
1815 // forward transform to fill the modal coeffs
1816 for (int i = 0; i < m_fields.size(); ++i)
1817 {
1818 m_fields[i]->SetPhysState(true);
1819 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1820 m_fields[i]->UpdateCoeffs());
1821 }
1822 }
1823 break;
1824
1825 case eTestSteadyZonal:
1826 {
1827 Array<OneD, NekDouble> eta0(nq);
1830 Array<OneD, NekDouble> zeta0(nq);
1831
1832 SteadyZonalFlow(0, eta0);
1833 m_fields[0]->SetPhys(eta0);
1834
1835 SteadyZonalFlow(1, u0);
1836 m_fields[1]->SetPhys(u0);
1837
1838 SteadyZonalFlow(2, v0);
1839 m_fields[2]->SetPhys(v0);
1840
1841 // ComputeVorticity(u0, v0, zeta0);
1842 m_Vorticity0 = m_fields[0]->Integral(zeta0);
1843
1844 m_Mass0 = ComputeMass(eta0);
1845 m_Energy0 = ComputeEnergy(eta0, u0, v0);
1846 m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1847
1848 // forward transform to fill the modal coeffs
1849 for (int i = 0; i < m_fields.size(); ++i)
1850 {
1851 m_fields[i]->SetPhysState(true);
1852 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1853 m_fields[i]->UpdateCoeffs());
1854 }
1855 }
1856 break;
1857
1858 case eTestUnsteadyZonal:
1859 {
1860 Array<OneD, NekDouble> eta0(nq);
1863
1864 UnsteadyZonalFlow(0, initialtime, eta0);
1865 m_fields[0]->SetPhys(eta0);
1866
1867 UnsteadyZonalFlow(1, initialtime, u0);
1868 m_fields[1]->SetPhys(u0);
1869
1870 UnsteadyZonalFlow(2, initialtime, v0);
1871 m_fields[2]->SetPhys(v0);
1872
1873 m_Mass0 = ComputeMass(eta0);
1874 m_Energy0 = ComputeEnergy(eta0, u0, v0);
1875 m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1876
1877 // forward transform to fill the modal coeffs
1878 for (int i = 0; i < m_fields.size(); ++i)
1879 {
1880 m_fields[i]->SetPhysState(true);
1881 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1882 m_fields[i]->UpdateCoeffs());
1883 }
1884 }
1885 break;
1886
1888 {
1889 Array<OneD, NekDouble> eta0(nq);
1892
1893 IsolatedMountainFlow(0, initialtime, eta0);
1894 m_fields[0]->SetPhys(eta0);
1895
1896 IsolatedMountainFlow(1, initialtime, u0);
1897 m_fields[1]->SetPhys(u0);
1898
1899 IsolatedMountainFlow(2, initialtime, v0);
1900 m_fields[2]->SetPhys(v0);
1901
1902 m_Mass0 = ComputeMass(eta0);
1903 m_Energy0 = ComputeEnergy(eta0, u0, v0);
1904 m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1905
1906 // forward transform to fill the modal coeffs
1907 for (int i = 0; i < m_fields.size(); ++i)
1908 {
1909 m_fields[i]->SetPhysState(true);
1910 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1911 m_fields[i]->UpdateCoeffs());
1912 }
1913 }
1914 break;
1915
1916 case eTestUnstableJet:
1917 {
1918 Array<OneD, NekDouble> eta0(nq);
1921
1922 UnstableJetFlow(0, initialtime, eta0);
1923 m_fields[0]->SetPhys(eta0);
1924
1925 UnstableJetFlow(1, initialtime, u0);
1926 m_fields[1]->SetPhys(u0);
1927
1928 UnstableJetFlow(2, initialtime, v0);
1929 m_fields[2]->SetPhys(v0);
1930
1931 m_Mass0 = ComputeMass(eta0);
1932 m_Energy0 = ComputeEnergy(eta0, u0, v0);
1933 m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1934
1935 // forward transform to fill the modal coeffs
1936 for (int i = 0; i < m_fields.size(); ++i)
1937 {
1938 m_fields[i]->SetPhysState(true);
1939 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1940 m_fields[i]->UpdateCoeffs());
1941 }
1942 }
1943 break;
1944
1945 case eTestRossbyWave:
1946 {
1947 Array<OneD, NekDouble> eta0(nq);
1950
1951 RossbyWave(0, eta0);
1952 m_fields[0]->SetPhys(eta0);
1953
1954 RossbyWave(1, u0);
1955 m_fields[1]->SetPhys(u0);
1956
1957 RossbyWave(2, v0);
1958 m_fields[2]->SetPhys(v0);
1959
1960 m_Mass0 = ComputeMass(eta0);
1961 m_Energy0 = ComputeEnergy(eta0, u0, v0);
1962 m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1963
1964 // forward transform to fill the modal coeffs
1965 for (int i = 0; i < m_fields.size(); ++i)
1966 {
1967 m_fields[i]->SetPhysState(true);
1968 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1969 m_fields[i]->UpdateCoeffs());
1970 }
1971 }
1972 break;
1973
1974 default:
1975 break;
1976 }
1977
1978 if (dumpInitialConditions)
1979 {
1980 // dump initial conditions to file
1981 std::string outname = m_sessionName + "_initial.chk";
1982 WriteFld(outname);
1983
1984 outname = m_sessionName + "_initialCART.chk";
1986 }
1987}
1988
1989void MMFSWE::TestSWE2Dproblem([[maybe_unused]] const NekDouble time,
1990 unsigned int field,
1991 Array<OneD, NekDouble> &outfield)
1992{
1993 int nq = m_fields[0]->GetNpoints();
1994
1998
1999 m_fields[0]->GetCoords(x0, x1, x2);
2000
2001 Array<OneD, NekDouble> eta0(nq);
2004
2006
2007 for (int i = 0; i < m_spacedim; ++i)
2008 {
2009 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2010 }
2011
2012 for (int i = 0; i < nq; ++i)
2013 {
2014 eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2015 (1.0 / cosh(0.395 * x0[i]))) *
2016 (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2017 exp(-0.5 * x1[i] * x1[i]);
2018 uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2019 (1.0 / cosh(0.395 * x0[i]))) *
2020 (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2021 exp(-0.5 * x1[i] * x1[i]);
2022 uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2023 (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2024 (1.0 / cosh(0.395 * x0[i]))) *
2025 (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2026 }
2027
2028 u0 = CartesianToMovingframes(uvec, 0);
2029 v0 = CartesianToMovingframes(uvec, 1);
2030
2031 switch (field)
2032 {
2033 case (0):
2034 {
2035 outfield = eta0;
2036 }
2037 break;
2038
2039 case (1):
2040 {
2041 outfield = u0;
2042 }
2043 break;
2044
2045 case (2):
2046 {
2047 outfield = v0;
2048 }
2049 break;
2050 }
2051}
2052
2054 Array<OneD, NekDouble> &outfield)
2055{
2056 int nq = GetTotPoints();
2057 NekDouble uhat, vhat;
2058 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2059 NekDouble x0j, x1j, x2j, tmp;
2060
2061 Array<OneD, NekDouble> eta(nq, 0.0);
2062 Array<OneD, NekDouble> u(nq, 0.0);
2063 Array<OneD, NekDouble> v(nq, 0.0);
2064
2066 for (int i = 0; i < m_spacedim; ++i)
2067 {
2068 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2069 }
2070
2074
2075 m_fields[0]->GetCoords(x, y, z);
2076
2077 for (int j = 0; j < nq; ++j)
2078 {
2079 x0j = x[j];
2080 x1j = y[j];
2081 x2j = z[j];
2082
2083 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2084 cos_theta);
2085
2086 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2087 // \sin \alpha + \sin \theta \cos \alpha )^2
2088 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
2089 sin_theta * cos(m_alpha);
2090 eta[j] = m_H0 - m_Hvar * tmp * tmp;
2091
2092 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2093 // || e^2 ||^2
2094 uhat = m_u0 * (cos_theta * cos(m_alpha) +
2095 sin_theta * cos_varphi * sin(m_alpha));
2096 vhat = -1.0 * m_u0 * sin_varphi * sin(m_alpha);
2097
2098 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2099 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2100 uvec[2][j] = vhat * cos_theta;
2101 }
2102
2103 // Projection of u onto the tangent plane with conserving the mag. of the
2104 // velocity.
2106
2107 for (int i = 0; i < m_spacedim; ++i)
2108 {
2109 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2110 }
2111
2112 // u is projected on the tangent plane with preserving its length
2113 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2114
2115 // Change it to the coordinate of moving frames
2116 // CartesianToMovingframes(0,uvecproj,u);
2117 // CartesianToMovingframes(1,uvecproj,v);
2118
2119 u = CartesianToMovingframes(uvec, 0);
2120 v = CartesianToMovingframes(uvec, 1);
2121
2122 switch (field)
2123 {
2124 case (0):
2125 {
2126 outfield = eta;
2127 }
2128 break;
2129
2130 case (1):
2131 {
2132 outfield = u;
2133 }
2134 break;
2135
2136 case (2):
2137 {
2138 outfield = v;
2139 }
2140 break;
2141 }
2142}
2143
2145{
2146 int nq = m_fields[0]->GetTotPoints();
2147
2148 Array<OneD, NekDouble> tmp(nq);
2149 Vmath::Vadd(nq, eta, 1, m_depth, 1, tmp, 1);
2150
2151 return m_fields[0]->Integral(tmp);
2152}
2153
2157{
2158 int nq = m_fields[0]->GetTotPoints();
2159
2160 Array<OneD, NekDouble> tmp(nq);
2161 Array<OneD, NekDouble> htmp(nq);
2162 Array<OneD, NekDouble> hstmp(nq);
2163
2164 Vmath::Vmul(nq, u, 1, u, 1, tmp, 1);
2165 Vmath::Vvtvp(nq, v, 1, v, 1, tmp, 1, tmp, 1);
2166 Vmath::Vmul(nq, eta, 1, tmp, 1, tmp, 1);
2167
2168 Vmath::Sadd(nq, m_H0, eta, 1, htmp, 1);
2169 Vmath::Vmul(nq, htmp, 1, htmp, 1, htmp, 1);
2170
2171 Vmath::Sadd(nq, -1.0 * m_H0, m_depth, 1, hstmp, 1);
2172 Vmath::Vmul(nq, hstmp, 1, hstmp, 1, hstmp, 1);
2173
2174 Vmath::Vsub(nq, htmp, 1, hstmp, 1, htmp, 1);
2175 Vmath::Smul(nq, m_g, htmp, 1, htmp, 1);
2176
2177 Vmath::Vadd(nq, htmp, 1, tmp, 1, tmp, 1);
2178 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2179
2180 return m_fields[0]->Integral(tmp);
2181}
2182
2186{
2187 int nq = m_fields[0]->GetTotPoints();
2188
2189 Array<OneD, NekDouble> hstartmp(nq);
2190 Array<OneD, NekDouble> tmp(nq);
2191
2192 Vmath::Vadd(nq, eta, 1, m_depth, 1, hstartmp, 1);
2193
2196 for (int i = 0; i < m_spacedim; ++i)
2197 {
2198 uvec[i] = Array<OneD, NekDouble>(nq);
2199 Curlu[i] = Array<OneD, NekDouble>(nq);
2200 }
2201
2202 ComputeVorticity(u, v, tmp);
2203
2204 Vmath::Vadd(nq, m_coriolis, 1, tmp, 1, tmp, 1);
2205 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
2206 Vmath::Vdiv(nq, tmp, 1, hstartmp, 1, tmp, 1);
2207 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2208
2209 return m_fields[0]->Integral(tmp);
2210}
2211
2212// Vorticity = \nabla v \cdot e^1 + v \nabla \cdot e^1 - ( \nabla u \cdot e^2 +
2213// u \nabla \cdot e^2 )
2216 Array<OneD, NekDouble> &Vorticity)
2217{
2218 int nq = m_fields[0]->GetTotPoints();
2219
2220 Array<OneD, NekDouble> tmp(nq);
2221
2222 Vorticity = Array<OneD, NekDouble>(nq, 0.0);
2223
2224 m_fields[0]->PhysDirectionalDeriv(m_movingframes[0], v, Vorticity);
2225 Vmath::Vvtvp(nq, &v[0], 1, &m_CurlMF[1][2][0], 1, &Vorticity[0], 1,
2226 &Vorticity[0], 1);
2227
2228 m_fields[0]->PhysDirectionalDeriv(m_movingframes[1], u, tmp);
2229 Vmath::Neg(nq, tmp, 1);
2230 Vmath::Vvtvp(nq, &u[0], 1, &m_CurlMF[0][2][0], 1, &tmp[0], 1, &tmp[0], 1);
2231
2232 Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2233}
2234
2236{
2237 int nq = m_fields[0]->GetNpoints();
2238
2239 Array<OneD, NekDouble> velcoeff(nq, 0.0);
2240
2241 Array<OneD, NekDouble> Dtmp0(nq);
2242 Array<OneD, NekDouble> Dtmp1(nq);
2243 Array<OneD, NekDouble> Dtmp2(nq);
2244 Array<OneD, NekDouble> Drv(nq);
2245
2246 vellc = Array<OneD, NekDouble>(nq, 0.0);
2247
2248 // m_vellc = \nabla m_vel \cdot tan_i
2249 Array<OneD, NekDouble> tmp(nq);
2250 Array<OneD, NekDouble> vessel(nq);
2251
2252 for (int j = 0; j < m_shapedim; ++j)
2253 {
2254 Vmath::Zero(nq, velcoeff, 1);
2255 for (int k = 0; k < m_spacedim; ++k)
2256 {
2257 // a_j = tan_j cdot m_vel
2258 Vmath::Vvtvp(nq, &m_movingframes[j][k * nq], 1, &m_velocity[k][0],
2259 1, &velcoeff[0], 1, &velcoeff[0], 1);
2260 }
2261
2262 // d a_j / d x^k
2263 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2264
2265 for (int k = 0; k < m_spacedim; ++k)
2266 {
2267 // tan_j^k ( d a_j / d x^k )
2268 switch (k)
2269 {
2270 case (0):
2271 {
2272 Vmath::Vvtvp(nq, &Dtmp0[0], 1, &m_movingframes[j][k * nq],
2273 1, &vellc[0], 1, &vellc[0], 1);
2274 }
2275 break;
2276
2277 case (1):
2278 {
2279 Vmath::Vvtvp(nq, &Dtmp1[0], 1, &m_movingframes[j][k * nq],
2280 1, &vellc[0], 1, &vellc[0], 1);
2281 }
2282 break;
2283
2284 case (2):
2285 {
2286 Vmath::Vvtvp(nq, &Dtmp2[0], 1, &m_movingframes[j][k * nq],
2287 1, &vellc[0], 1, &vellc[0], 1);
2288 }
2289 break;
2290 }
2291 }
2292 }
2293}
2294
2295void MMFSWE::UnsteadyZonalFlow(unsigned int field, const NekDouble time,
2296 Array<OneD, NekDouble> &outfield)
2297{
2298 int nq = GetTotPoints();
2299 NekDouble uhat, vhat;
2300 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2301 NekDouble x0j, x1j, x2j, tmp;
2302
2303 NekDouble TR, Ttheta;
2304
2305 Array<OneD, NekDouble> eta(nq, 0.0);
2306 Array<OneD, NekDouble> u(nq, 0.0);
2307 Array<OneD, NekDouble> v(nq, 0.0);
2308
2310 for (int i = 0; i < m_spacedim; ++i)
2311 {
2312 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2313 }
2314
2318
2319 m_fields[0]->GetCoords(x, y, z);
2320
2321 for (int j = 0; j < nq; ++j)
2322 {
2323 x0j = x[j];
2324 x1j = y[j];
2325 x2j = z[j];
2326
2327 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2328 cos_theta);
2329
2330 // \eta = ( - ( u_0 ( - T_R sin \alpha \cos \theta + \cos \alpha \sin
2331 // \theta ) + \Omega \sin \theta )^2 + \Omega \sin \theta )^2 + (\Omega
2332 // \sin \theta )^2
2333 TR =
2334 cos_varphi * cos(m_Omega * time) - sin_varphi * sin(m_Omega * time);
2335 Ttheta =
2336 sin_varphi * cos(m_Omega * time) + cos_varphi * sin(m_Omega * time);
2337 tmp = -1.0 * TR * sin(m_alpha) * cos_theta + cos(m_alpha) * sin_theta;
2338
2339 eta[j] = -1.0 * (m_u0 * tmp + m_Omega * sin_theta) *
2340 (m_u0 * tmp + m_Omega * sin_theta) +
2341 m_Omega * m_Omega * sin_theta * sin_theta;
2342 eta[j] = 0.5 * eta[j] / m_g;
2343
2344 // u = u_0*(TR*\sin \alpha * \sin \theta + \cos \alpha * \cos \theta
2345 // v = - u_0 Ttheta * \sin \alpha
2346 uhat =
2347 m_u0 * (TR * sin(m_alpha) * sin_theta + cos(m_alpha) * cos_theta);
2348 vhat = -1.0 * m_u0 * Ttheta * sin(m_alpha);
2349
2350 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2351 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2352 uvec[2][j] = vhat * cos_theta;
2353 }
2354
2355 // Projection of u onto the tangent plane with conserving the mag. of the
2356 // velocity.
2358
2359 for (int i = 0; i < m_spacedim; ++i)
2360 {
2361 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2362 }
2363
2364 // u is projected on the tangent plane with preserving its length
2365 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2366
2367 // Change it to the coordinate of moving frames
2368 // CartesianToMovingframes(0,uvecproj,u);
2369 // CartesianToMovingframes(1,uvecproj,v);
2370
2371 u = CartesianToMovingframes(uvec, 0);
2372 v = CartesianToMovingframes(uvec, 1);
2373
2374 switch (field)
2375 {
2376 case (0):
2377 {
2378 outfield = eta;
2379 }
2380 break;
2381
2382 case (1):
2383 {
2384 outfield = u;
2385 }
2386 break;
2387
2388 case (2):
2389 {
2390 outfield = v;
2391 }
2392 break;
2393 }
2394}
2395
2397 [[maybe_unused]] const NekDouble time,
2398 Array<OneD, NekDouble> &outfield)
2399{
2400 int nq = GetTotPoints();
2401
2402 NekDouble uhat, vhat;
2403 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2404 NekDouble x0j, x1j, x2j;
2405
2406 Array<OneD, NekDouble> eta(nq, 0.0);
2407 Array<OneD, NekDouble> u(nq, 0.0);
2408 Array<OneD, NekDouble> v(nq, 0.0);
2409
2411 for (int i = 0; i < m_spacedim; ++i)
2412 {
2413 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2414 }
2415
2419
2420 m_fields[0]->GetCoords(x, y, z);
2421
2422 for (int j = 0; j < nq; ++j)
2423 {
2424 x0j = x[j];
2425 x1j = y[j];
2426 x2j = z[j];
2427
2428 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2429 cos_theta);
2430
2431 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2432 eta[j] = (-1.0 / m_g) * (m_Omega * m_u0 + 0.5 * m_u0 * m_u0) *
2433 sin_theta * sin_theta;
2434
2435 uhat = m_u0 * cos_theta;
2436 vhat = 0.0;
2437
2438 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2439 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2440 uvec[2][j] = vhat * cos_theta;
2441 }
2442
2443 // Projection of u onto the tangent plane with conserving the mag. of the
2444 // velocity.
2446
2447 for (int i = 0; i < m_spacedim; ++i)
2448 {
2449 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2450 }
2451
2452 // u is projected on the tangent plane with preserving its length
2453 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2454
2455 // Change it to the coordinate of moving frames
2456 // CartesianToMovingframes(0,uvecproj,u);
2457 // CartesianToMovingframes(1,uvecproj,v);
2458
2459 u = CartesianToMovingframes(uvec, 0);
2460 v = CartesianToMovingframes(uvec, 1);
2461
2462 switch (field)
2463 {
2464 case (0):
2465 {
2466 outfield = eta;
2467 }
2468 break;
2469
2470 case (1):
2471 {
2472 outfield = u;
2473 }
2474 break;
2475
2476 case (2):
2477 {
2478 outfield = v;
2479 }
2480 break;
2481 }
2482}
2483
2485 [[maybe_unused]] const NekDouble time,
2486 Array<OneD, NekDouble> &outfield)
2487{
2488 int nq = GetTotPoints();
2489
2490 NekDouble uhat, vhat;
2491 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2492 NekDouble x0j, x1j, x2j;
2493 NekDouble Ttheta, Tphi;
2494
2495 Array<OneD, NekDouble> eta(nq, 0.0);
2496 Array<OneD, NekDouble> u(nq, 0.0);
2497 Array<OneD, NekDouble> v(nq, 0.0);
2498
2500 for (int i = 0; i < m_spacedim; ++i)
2501 {
2502 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2503 }
2504
2508
2509 m_fields[0]->GetCoords(x, y, z);
2510
2511 int Nint = 1000;
2512 NekDouble dth, intj;
2513 for (int j = 0; j < nq; ++j)
2514 {
2515 x0j = x[j];
2516 x1j = y[j];
2517 x2j = z[j];
2518
2519 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2520 cos_theta);
2521
2522 Ttheta = atan2(sin_theta, cos_theta);
2523 Tphi = atan2(sin_varphi, cos_varphi);
2524
2525 uhat = ComputeUnstableJetuphi(Ttheta);
2526 vhat = 0.0;
2527
2528 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2529 dth = Ttheta / Nint;
2530 eta[j] = dth * 0.5 *
2532 for (int i = 1; i < Nint - 1; i++)
2533 {
2534 intj = i * dth;
2535 eta[j] = eta[j] + dth * ComputeUnstableJetEta(intj);
2536 }
2537 eta[j] = (-1.0 / m_g) * eta[j];
2538
2539 // Add perturbation
2540 if (m_PurturbedJet)
2541 {
2542 eta[j] = eta[j] + m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2543 exp(-225.0 * (m_pi / 4.0 - Ttheta) *
2544 (m_pi / 4.0 - Ttheta));
2545 }
2546
2547 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2548 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2549 uvec[2][j] = vhat * cos_theta;
2550 }
2551
2552 // Projection of u onto the tangent plane with conserving the mag. of the
2553 // velocity.
2555
2556 for (int i = 0; i < m_spacedim; ++i)
2557 {
2558 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2559 }
2560
2561 // u is projected on the tangent plane with preserving its length
2562 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2563
2564 // Change it to the coordinate of moving frames
2565 // CartesianToMovingframes(0,uvecproj,u);
2566 // CartesianToMovingframes(1,uvecproj,v);
2567
2568 u = CartesianToMovingframes(uvec, 0);
2569 v = CartesianToMovingframes(uvec, 1);
2570
2571 switch (field)
2572 {
2573 case (0):
2574 {
2575 outfield = eta;
2576 }
2577 break;
2578
2579 case (1):
2580 {
2581 outfield = u;
2582 }
2583 break;
2584
2585 case (2):
2586 {
2587 outfield = v;
2588 }
2589 break;
2590 }
2591}
2592
2594{
2595 int nq = GetTotPoints();
2596 NekDouble uhat, vhat;
2597 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2598 NekDouble x0j, x1j, x2j;
2599 NekDouble Ath, Bth, Cth, tmp;
2600
2601 Array<OneD, NekDouble> eta(nq, 0.0);
2602 Array<OneD, NekDouble> u(nq, 0.0);
2603 Array<OneD, NekDouble> v(nq, 0.0);
2604
2606 for (int i = 0; i < m_spacedim; ++i)
2607 {
2608 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2609 }
2610
2614
2615 m_fields[0]->GetCoords(x, y, z);
2616
2617 NekDouble R = 4.0;
2618 NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2619 NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2620
2621 // disturbancees of Rossby-Haurwitz Wave
2622 NekDouble x0d, y0d, z0d, phi0, theta0;
2623 // NekDouble rad_earth = 6.37122 * 1000000;
2624
2625 phi0 = 40.0 * m_pi / 180.0;
2626 theta0 = 50.0 * m_pi / 180.0;
2627
2628 x0d = cos(phi0) * cos(theta0);
2629 y0d = sin(phi0) * cos(theta0);
2630 z0d = sin(theta0);
2631
2632 for (int j = 0; j < nq; ++j)
2633 {
2634 x0j = x[j];
2635 x1j = y[j];
2636 x2j = z[j];
2637
2638 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2639 cos_theta);
2640
2641 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2642 // \sin \alpha + \sin \theta \cos \alpha )^2
2643
2644 // tmp = cos^{2R} \theta, R = 4;
2645 cos2theta = cos_theta * cos_theta;
2646 cosRm1theta = cos_theta * cos2theta;
2647 cosRtheta = cos2theta * cos2theta;
2648 cos6theta = cos2theta * cosRtheta;
2649 cos2Rtheta = cosRtheta * cosRtheta;
2650
2651 tmp = (0.5 * m_angfreq) * (2.0 * m_Omega + m_angfreq);
2652 Ath = tmp * cos2theta +
2653 0.25 * m_K * m_K * cos6theta *
2654 ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2655 2 * R * R);
2656
2657 tmp = (2.0 * m_K) * (m_Omega + m_angfreq) / ((R + 1.0) * (R + 2.0));
2658 Bth = tmp * cosRtheta *
2659 ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2660
2661 Cth =
2662 0.25 * m_K * m_K * cos2Rtheta * ((R + 1.0) * cos2theta - (R + 2.0));
2663
2664 // cos (2 \phi) = 2 * cos \phi * cos \phi - 1.0
2665 cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2666 cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2667 cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2668
2669 // sin (2 \phi) = 2 * cos \phi * sin \phi
2670 sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2671
2672 eta[j] = m_H0 + (1.0 / m_g) * (Ath + Bth * cos4phi + Cth * cos8phi);
2673
2674 // disturbances is added
2676 {
2677 eta[j] = eta[j] *
2678 (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2679 }
2680
2681 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2682 // || e^2 ||^2
2683 uhat = m_angfreq * cos_theta +
2684 m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2685 cos4phi;
2686 vhat = -1.0 * m_K * R * cosRm1theta * sin_theta * sin4phi;
2687
2688 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2689 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2690 uvec[2][j] = vhat * cos_theta;
2691 }
2692
2693 // NekDouble etamin, etaaver;
2694 // etamin = Vmath::Vmin(nq, eta, 1)*rad_earth;
2695 // etaaver = Vmath::Vsum(nq, eta, 1)/nq*rad_earth;
2696
2697 // Projection of u onto the tangent plane with conserving the mag. of the
2698 // velocity.
2700
2701 for (int i = 0; i < m_spacedim; ++i)
2702 {
2703 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2704 }
2705
2706 // u is projected on the tangent plane with preserving its length
2707 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2708
2709 // Change it to the coordinate of moving frames
2710 // CartesianToMovingframes(0,uvecproj,u);
2711 // CartesianToMovingframes(1,uvecproj,v);
2712
2713 u = CartesianToMovingframes(uvec, 0);
2714 v = CartesianToMovingframes(uvec, 1);
2715
2716 switch (field)
2717 {
2718 case (0):
2719 {
2720 outfield = eta;
2721 }
2722 break;
2723
2724 case (1):
2725 {
2726 outfield = u;
2727 }
2728 break;
2729
2730 case (2):
2731 {
2732 outfield = v;
2733 }
2734 break;
2735 }
2736}
2737
2739{
2740 NekDouble uphi, f, dh;
2741
2742 uphi = ComputeUnstableJetuphi(theta);
2743 f = 2.0 * m_Omega * sin(theta);
2744
2745 dh = f * uphi + tan(theta) * uphi * uphi;
2746
2747 return dh;
2748}
2749
2751{
2752 NekDouble uphi;
2753
2754 if ((theta > m_theta0) && (theta < m_theta1))
2755 {
2756 uphi = (m_uthetamax / m_en) *
2757 exp(1.0 / (theta - m_theta0) / (theta - m_theta1));
2758 }
2759
2760 else
2761 {
2762 uphi = 0.0;
2763 }
2764
2765 return uphi;
2766}
2767
2769{
2770 int nq = m_fields[0]->GetTotPoints();
2771 int ncoeffs = m_fields[0]->GetNcoeffs();
2772
2773 NekDouble rad_earth = 6.37122 * 1000000;
2774
2775 int nvariables = 7;
2776
2777 // vector u in Cartesian coordinates
2778 std::vector<std::string> variables(nvariables);
2779
2780 variables[0] = "eta";
2781 variables[1] = "hstar";
2782 variables[2] = "vorticity";
2783 variables[3] = "ux";
2784 variables[4] = "uy";
2785 variables[5] = "uz";
2786 variables[6] = "null";
2787
2788 // Obtain \vec{u} in cartesian coordinate
2789 Array<OneD, Array<OneD, NekDouble>> fieldphys(nvariables);
2790 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2791 for (int i = 0; i < nvariables; ++i)
2792 {
2793 fieldphys[i] = Array<OneD, NekDouble>(nq, 0.0);
2794 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2795 }
2796
2797 Vmath::Smul(nq, rad_earth, &(m_fields[0]->GetPhys())[0], 1,
2798 &fieldphys[0][0], 1);
2799 // Vmath::Vadd(nq, &(m_fields[0]->GetPhys())[0], 1, &m_depth[0], 1,
2800 // &fieldphys[1][0], 1);
2801
2802 Vmath::Vcopy(nq, &m_depth[0], 1, &fieldphys[1][0], 1);
2803 Vmath::Neg(nq, &fieldphys[1][0], 1);
2804 Vmath::Sadd(nq, m_H0, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2805 Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2806
2807 Array<OneD, NekDouble> utmp(nq);
2808 Array<OneD, NekDouble> vtmp(nq);
2809
2810 Vmath::Vcopy(nq, &(m_fields[1]->GetPhys())[0], 1, &utmp[0], 1);
2811 Vmath::Vcopy(nq, &(m_fields[2]->GetPhys())[0], 1, &vtmp[0], 1);
2812
2813 ComputeVorticity(utmp, vtmp, fieldphys[2]);
2814 // u_x = u e^1_x + v e^2_x
2815 int indx = 3;
2816 for (int k = 0; k < m_spacedim; ++k)
2817 {
2818 Vmath::Vmul(nq, &utmp[0], 1, &m_movingframes[0][k * nq], 1,
2819 &fieldphys[k + indx][0], 1);
2820 Vmath::Vvtvp(nq, &vtmp[0], 1, &m_movingframes[1][k * nq], 1,
2821 &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2822 }
2823
2824 for (int i = 0; i < nvariables; ++i)
2825 {
2826 m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2827 }
2828
2829 WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2830}
2831
2832// (h, hu, hv) -> (\eta, u, v)
2834 const Array<OneD, const Array<OneD, NekDouble>> &physin,
2836{
2837 int nq = GetTotPoints();
2838
2839 if (physin[0].data() == physout[0].data())
2840 {
2841 // copy indata and work with tmp array
2843 for (int i = 0; i < 3; ++i)
2844 {
2845 // deep copy
2846 tmp[i] = Array<OneD, NekDouble>(nq);
2847 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2848 }
2849
2850 // \eta = h - d
2851 Vmath::Vsub(nq, tmp[0], 1, m_depth, 1, physout[0], 1);
2852
2853 // u = hu/h
2854 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2855
2856 // v = hv/h
2857 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2858 }
2859 else
2860 {
2861 // \eta = h - d
2862 Vmath::Vsub(nq, physin[0], 1, m_depth, 1, physout[0], 1);
2863
2864 // u = hu/h
2865 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2866
2867 // v = hv/h
2868 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2869 }
2870}
2871
2873 const Array<OneD, const Array<OneD, NekDouble>> &physin,
2875{
2876
2877 int nq = GetTotPoints();
2878
2879 if (physin[0].data() == physout[0].data())
2880 {
2881 // copy indata and work with tmp array
2883 for (int i = 0; i < 3; ++i)
2884 {
2885 // deep copy
2886 tmp[i] = Array<OneD, NekDouble>(nq);
2887 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2888 }
2889
2890 // h = \eta + d
2891 Vmath::Vadd(nq, tmp[0], 1, m_depth, 1, physout[0], 1);
2892
2893 // hu = h * u
2894 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2895
2896 // hv = h * v
2897 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2898 }
2899 else
2900 {
2901 // h = \eta + d
2902 Vmath::Vadd(nq, physin[0], 1, m_depth, 1, physout[0], 1);
2903
2904 // hu = h * u
2905 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2906
2907 // hv = h * v
2908 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2909 }
2910}
2911
2912// To obtain [eta, u, v ] from [H, Hu, Hv]
2914{
2915 int nq = GetTotPoints();
2916
2917 // u = hu/h
2918 Vmath::Vdiv(nq, m_fields[1]->GetPhys(), 1, m_fields[0]->GetPhys(), 1,
2919 m_fields[1]->UpdatePhys(), 1);
2920
2921 // v = hv/ v
2922 Vmath::Vdiv(nq, m_fields[2]->GetPhys(), 1, m_fields[0]->GetPhys(), 1,
2923 m_fields[2]->UpdatePhys(), 1);
2924
2925 // \eta = h - d
2926 Vmath::Vsub(nq, m_fields[0]->GetPhys(), 1, m_depth, 1,
2927 m_fields[0]->UpdatePhys(), 1);
2928}
2929
2931{
2932 int nq = GetTotPoints();
2933
2934 // h = \eta + d
2935 Vmath::Vadd(nq, m_fields[0]->GetPhys(), 1, m_depth, 1,
2936 m_fields[0]->UpdatePhys(), 1);
2937
2938 // hu = h * u
2939 Vmath::Vmul(nq, m_fields[0]->GetPhys(), 1, m_fields[1]->GetPhys(), 1,
2940 m_fields[1]->UpdatePhys(), 1);
2941
2942 // hv = h * v
2943 Vmath::Vmul(nq, m_fields[0]->GetPhys(), 1, m_fields[2]->GetPhys(), 1,
2944 m_fields[2]->UpdatePhys(), 1);
2945}
2946
2948{
2949 // Construct beta
2950 int i, k;
2951 int n = 1, m = 1;
2952 int nq = m_fields[0]->GetTotPoints();
2953
2954 NekDouble alpha, beta_theta, beta_phi;
2955
2956 NekDouble xp, yp, zp, Re;
2957 NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2958 NekDouble cosntheta3;
2959
2960 NekDouble thetax, thetay, thetaz;
2961 NekDouble phix, phiy, phiz;
2962
2966
2969
2970 Array<OneD, NekDouble> vorticitycompt(nq);
2971 Array<OneD, NekDouble> vorticityexact(nq);
2972
2973 m_fields[0]->GetCoords(x, y, z);
2974
2976 for (i = 0; i < m_spacedim; ++i)
2977 {
2978 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2979 }
2980
2981 // Generate SphericalCoords
2982 for (k = 0; k < nq; ++k)
2983 {
2984 xp = x[k];
2985 yp = y[k];
2986 zp = z[k];
2987
2988 Re = sqrt(xp * xp + yp * yp + zp * zp);
2989
2990 CartesianToSpherical(xp, yp, zp, sin_varphi, cos_varphi, sin_theta,
2991 cos_theta);
2992
2993 alpha = sin_varphi;
2994
2995 theta = atan2(sin_theta, cos_theta);
2996 phi = atan2(sin_varphi, cos_varphi);
2997
2998 cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
2999
3000 beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
3001 beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3002
3003 thetax = -1.0 * cos_varphi * sin_theta;
3004 thetay = -1.0 * sin_varphi * sin_theta;
3005 thetaz = cos_theta;
3006
3007 phix = -1.0 * sin_varphi;
3008 phiy = cos_varphi;
3009 phiz = 0.0;
3010
3011 uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3012 uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3013 uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3014
3015 vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3016 cos_varphi * cos(m * phi) * sin(n * theta);
3017 }
3018
3019 u = CartesianToMovingframes(uvec, 0);
3020 v = CartesianToMovingframes(uvec, 1);
3021
3022 std::cout << "chi migi1" << std::endl;
3023
3024 ComputeVorticity(u, v, vorticitycompt);
3025
3026 Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3027
3028 std::cout << "Vorticity: L2 error = " << AvgAbsInt(vorticitycompt)
3029 << ", Linf error = " << Vmath::Vamax(nq, vorticitycompt, 1)
3030 << std::endl;
3031}
3032
3034 unsigned int field,
3035 [[maybe_unused]] const Array<OneD, NekDouble> &exactsoln, bool Normalised)
3036{
3037 int nq = m_fields[field]->GetNpoints();
3038 NekDouble L2error = -1.0;
3039
3040 if (m_NumQuadPointsError == 0)
3041 {
3042 if (m_fields[field]->GetPhysState() == false)
3043 {
3044 m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
3045 m_fields[field]->UpdatePhys());
3046 }
3047
3048 switch (field)
3049 {
3050 case (0):
3051 {
3052 // I(h) = I (h - h_exact ) / I (h_exact)
3053 Array<OneD, NekDouble> exactsolution(nq);
3054 v_EvaluateExactSolution(0, exactsolution, m_time);
3055
3056 // exactsoln = u - u_T so that L2 compute u_T
3057 NekDouble L2exact = m_fields[0]->Integral(exactsolution);
3058
3059 Vmath::Vsub(nq, &(m_fields[0]->GetPhys())[0], 1,
3060 &exactsolution[0], 1, &exactsolution[0], 1);
3061 Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3062
3063 L2error = (m_fields[0]->Integral(exactsolution)) / L2exact;
3064 }
3065 break;
3066
3067 case (1):
3068 {
3069 // I2 (u) = I( (u - u_ext)^2 + (v - v_ext)^2 )^{1/2} / I(
3070 // u_ext^2 + v_ext^2 )^{1/2}
3071 Array<OneD, NekDouble> exactu(nq);
3072 Array<OneD, NekDouble> exactv(nq);
3073 Array<OneD, NekDouble> tmp(nq);
3074
3075 // L2exact = \int (\sqrt{exactu*exactu+exactv*exactv})
3076 NekDouble L2exact;
3077 v_EvaluateExactSolution(1, exactu, m_time);
3078 v_EvaluateExactSolution(2, exactv, m_time);
3079 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3080 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3081 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3082
3083 L2exact = m_fields[1]->Integral(tmp);
3084
3085 // L2exact = \int
3086 // (\sqrt{(u-exactu)*(u-exactu)+(v-exactv)*(v-exactv)})
3087 Vmath::Vsub(nq, &(m_fields[1]->GetPhys())[0], 1, &exactu[0], 1,
3088 &exactu[0], 1);
3089 Vmath::Vsub(nq, &(m_fields[2]->GetPhys())[0], 1, &exactv[0], 1,
3090 &exactv[0], 1);
3091 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3092 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3093 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3094
3095 L2error = (m_fields[1]->Integral(tmp)) / L2exact;
3096 }
3097 break;
3098
3099 case (2):
3100 {
3101 L2error = 0.0;
3102 }
3103 break;
3104
3105 default:
3106 break;
3107 }
3108
3109 if (Normalised == true)
3110 {
3112
3113 NekDouble Vol = m_fields[field]->Integral(one);
3114 m_comm->AllReduce(Vol, LibUtilities::ReduceSum);
3115
3116 L2error = sqrt(L2error * L2error / Vol);
3117 }
3118 }
3119 else
3120 {
3121 Array<OneD, NekDouble> L2INF(2);
3122 L2INF = ErrorExtraPoints(field);
3123 L2error = L2INF[0];
3124 }
3125
3126 return L2error;
3127}
3128
3130 unsigned int field,
3131 [[maybe_unused]] const Array<OneD, NekDouble> &exactsoln)
3132{
3133 NekDouble LinfError = -1.0;
3134
3135 if (m_fields[field]->GetPhysState() == false)
3136 {
3137 m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
3138 m_fields[field]->UpdatePhys());
3139 }
3140
3141 int nq = m_fields[field]->GetNpoints();
3142
3143 // Obtain \vec{u} in cartesian coordinate
3147
3148 m_fields[0]->GetCoords(x, y, z);
3149
3150 switch (field)
3151 {
3152 case (0):
3153 {
3154 NekDouble Letaint;
3155
3156 Array<OneD, NekDouble> exactsolution(nq);
3157
3158 EvaluateExactSolution(field, exactsolution, m_time);
3159 LinfError = m_fields[field]->Linf(m_fields[field]->GetPhys(),
3160 exactsolution);
3161
3162 Letaint = Vmath::Vamax(nq, exactsolution, 1);
3163
3164 Vmath::Vsub(nq, &(m_fields[0]->GetPhys())[0], 1, &exactsolution[0],
3165 1, &exactsolution[0], 1);
3166
3167 LinfError = fabs(LinfError / Letaint);
3168 }
3169 break;
3170
3171 case (1):
3172 {
3173 Array<OneD, NekDouble> exactu(nq);
3174 Array<OneD, NekDouble> exactv(nq);
3175 Array<OneD, NekDouble> tmpu(nq);
3176 Array<OneD, NekDouble> tmpv(nq);
3177 Array<OneD, NekDouble> Lerr(nq);
3179
3180 EvaluateExactSolution(1, exactu, m_time);
3181 EvaluateExactSolution(2, exactv, m_time);
3182
3183 // Compute max[sqrt{(u-uex)^2 + (v-vex)^2}]
3184 Vmath::Vcopy(nq, &(m_fields[1]->UpdatePhys())[0], 1, &tmpu[0], 1);
3185 Vmath::Vcopy(nq, &(m_fields[2]->UpdatePhys())[0], 1, &tmpv[0], 1);
3186
3187 Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3188 Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3189
3190 Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3191 Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3192
3193 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3194 Vmath::Vsqrt(nq, &Lerr[0], 1, &Lerr[0], 1);
3195
3196 // uT = max[sqrt( u_T^2 + v_T^2 ) ]
3197 Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3198 Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3199 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &uT[0], 1);
3200 Vmath::Vsqrt(nq, &uT[0], 1, &uT[0], 1);
3201
3202 LinfError = Vmath::Vamax(nq, Lerr, 1) / Vmath::Vamax(nq, uT, 1);
3203 }
3204 break;
3205
3206 case (2):
3207 {
3208 LinfError = 0.0;
3209 }
3210 break;
3211
3212 default:
3213 break;
3214 }
3215
3216 return LinfError;
3217}
3218
3220 Array<OneD, NekDouble> &outfield,
3221 const NekDouble time)
3222{
3223 switch (m_TestType)
3224 {
3225 case eTestSteadyZonal:
3226 {
3227 SteadyZonalFlow(field, outfield);
3228 }
3229 break;
3230
3231 case eTestUnsteadyZonal:
3232 {
3233 UnsteadyZonalFlow(field, time, outfield);
3234 }
3235 break;
3236
3238 {
3239 IsolatedMountainFlow(field, time, outfield);
3240 }
3241 break;
3242
3243 case eTestUnstableJet:
3244 {
3245 UnstableJetFlow(field, time, outfield);
3246 }
3247 break;
3248
3249 case eTestRossbyWave:
3250 {
3251 RossbyWave(field, outfield);
3252 }
3253 break;
3254
3255 case eTestPlane:
3256 {
3257 TestSWE2Dproblem(time, field, outfield);
3258 }
3259 break;
3260
3261 default:
3262 break;
3263 }
3264}
3265
3267{
3270
3271 switch (m_TestType)
3272 {
3273 case eTestSteadyZonal:
3274 {
3275 SolverUtils::AddSummaryItem(s, "Rotation Angle", m_alpha);
3276 }
3277 break;
3278
3279 case eTestRossbyWave:
3280 {
3281 SolverUtils::AddSummaryItem(s, "RossbyDistrubance",
3283 }
3284 break;
3285
3286 case eTestUnstableJet:
3287 {
3288 SolverUtils::AddSummaryItem(s, "PurturbedJet", m_PurturbedJet);
3289 }
3290 break;
3291
3292 default:
3293 break;
3294 }
3295}
3296
3297} // end namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
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 ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2214
void AddDivForGradient(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:543
void AddCoriolis(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1207
void UnstableJetFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2484
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
Definition: MMFSWE.cpp:62
int m_PurturbedJet
Definition: MMFSWE.h:268
NekDouble m_k2
Definition: MMFSWE.h:92
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
Definition: MMFSWE.cpp:2750
MMFSWE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: MMFSWE.cpp:52
void v_SetInitialConditions(const NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0) override
Definition: MMFSWE.cpp:1792
NekDouble m_alpha
Definition: MMFSWE.h:91
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFSWE.h:66
void WeakDGSWEDirDeriv(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField)
Definition: MMFSWE.cpp:473
NekDouble m_en
Definition: MMFSWE.h:94
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2593
void ComputeNablaCdotVelocity(Array< OneD, NekDouble > &vellc)
Definition: MMFSWE.cpp:2235
NekDouble v_L2Error(unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised) override
Virtual function for the L_2 error computation between fields and a given exact solution.
Definition: MMFSWE.cpp:3033
NekDouble m_Energy0
Definition: MMFSWE.h:97
void EvaluateStandardCoriolis(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1765
static std::string className
Name of class.
Definition: MMFSWE.h:76
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSWE.cpp:573
NekDouble m_theta0
Definition: MMFSWE.h:94
void PrimitiveToConservative()
Definition: MMFSWE.cpp:2930
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Definition: MMFSWE.cpp:3219
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &inarray, NekDouble time)
Definition: MMFSWE.cpp:1383
void IsolatedMountainFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2396
NekDouble m_hbar
Definition: MMFSWE.h:94
NekDouble m_K
Definition: MMFSWE.h:95
Array< OneD, NekDouble > m_coriolis
Coriolis force.
Definition: MMFSWE.h:86
void ComputeMagAndDot(const int index, NekDouble &MageF1, NekDouble &MageF2, NekDouble &MageB1, NekDouble &MageB2, NekDouble &eF1_cdot_eB1, NekDouble &eF1_cdot_eB2, NekDouble &eF2_cdot_eB1, NekDouble &eF2_cdot_eB2)
Definition: MMFSWE.cpp:1169
NekDouble m_Hvar
Definition: MMFSWE.h:91
void AddRotation(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1281
TestType m_TestType
Definition: MMFSWE.h:79
void AverageFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:962
NekDouble m_Vorticity0
Definition: MMFSWE.h:97
void AddElevationEffect(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1252
void Checkpoint_Output_Cartesian(std::string outname)
Definition: MMFSWE.cpp:2768
void Computehhuhvflux(NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, NekDouble hstar, NekDouble &hflux, NekDouble &huflux, NekDouble &hvflux)
Definition: MMFSWE.cpp:884
NekDouble v_LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln) override
Virtual function for the L_inf error computation between fields and a given exact solution.
Definition: MMFSWE.cpp:3129
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
Definition: MMFSWE.cpp:419
void TestVorticityComputation()
Definition: MMFSWE.cpp:2947
void NumericalSWEFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
Definition: MMFSWE.cpp:668
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
Definition: MMFSWE.cpp:1364
void RiemannSolverHLLC(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:839
NekDouble ComputeEnstrophy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2183
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
Definition: MMFSWE.cpp:3266
NekDouble m_uthetamax
Definition: MMFSWE.h:94
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
Definition: MMFSWE.cpp:2144
void EvaluateWaterDepth(void)
Definition: MMFSWE.cpp:1563
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2295
NekDouble ComputeEnergy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2154
void LaxFriedrichFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:1016
NekDouble ComputeUnstableJetEta(const NekDouble theta)
Definition: MMFSWE.cpp:2738
NekDouble m_Mass0
Definition: MMFSWE.h:97
void EvaluateCoriolis(void)
Definition: MMFSWE.cpp:1703
void EvaluateCoriolisForZonalFlow(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1733
void v_DoInitialise(bool dumpInitialConditions=false) override
Sets up initial conditions.
Definition: MMFSWE.cpp:1546
NekDouble m_theta1
Definition: MMFSWE.h:94
int m_AddRotation
Definition: MMFSWE.h:88
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
Definition: MMFSWE.h:103
int m_planeNumber
Definition: MMFSWE.h:109
void v_DoSolve() override
Solves an unsteady problem.
Definition: MMFSWE.cpp:242
NekDouble m_angfreq
Definition: MMFSWE.h:95
void Compute_demdt_cdot_ek(const int indm, const int indk, const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1328
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Definition: MMFSWE.cpp:1433
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:1989
NekDouble m_Omega
Definition: MMFSWE.h:91
int m_RossbyDisturbance
Definition: MMFSWE.h:267
NekDouble m_H0
Definition: MMFSWE.h:91
int m_AddCoriolis
Definition: MMFSWE.h:88
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
Definition: MMFSWE.h:83
NekDouble m_u0
Definition: MMFSWE.h:91
void ConservativeToPrimitive()
Definition: MMFSWE.cpp:2913
NekDouble m_Enstrophy0
Definition: MMFSWE.h:97
void SteadyZonalFlow(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2053
void RusanovFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:1091
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:82
NekDouble m_g
Definition: MMFSWE.h:90
NekDouble m_hs0
Definition: MMFSWE.h:93
int m_spacedim
Spatial dimension (>= expansion dim).
int m_expdim
Expansion dimension.
LibUtilities::CommSharedPtr m_comm
Communicator.
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.
SOLVER_UTILS_EXPORT int GetExpSize()
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints(unsigned int field)
Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf].
SOLVER_UTILS_EXPORT int GetTotPoints()
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
int m_steps
Number of steps to take.
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SOLVER_UTILS_EXPORT NekDouble LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Linf error computation.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
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
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:192
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:189
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Virtual function for generating summary information.
Definition: MMFSystem.cpp:2463
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
Definition: MMFSystem.h:190
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2307
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:183
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Definition: MMFSystem.h:196
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
Definition: MMFSystem.h:197
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble > > &uvec, unsigned int field)
Definition: MMFSystem.cpp:771
SOLVER_UTILS_EXPORT void CopyBoundaryTrace(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
Definition: MMFSystem.cpp:835
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:187
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
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
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:193
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.
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:135
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekZeroTol
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
EquationSystemFactory & GetEquationSystemFactory()
@ eAverage
averaged (or centred) flux
Definition: MMFSystem.h:77
@ eLaxFriedrich
Lax-Friedrich flux.
Definition: MMFSystem.h:78
@ eHLLC
Harten-Lax-Leer Contact wave flux.
Definition: MMFSystem.h:82
@ eRusanov
Upwind.
Definition: MMFSystem.h:80
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
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::vector< double > z(NPUPPER)
@ eTestUnstableJet
Definition: MMFSWE.h:49
@ eTestSteadyZonal
Definition: MMFSWE.h:46
@ eTestPlane
Definition: MMFDiffusion.h:48
@ eTestUnsteadyZonal
Definition: MMFSWE.h:47
@ SIZE_TestType
Length of enum list.
Definition: MMFDiffusion.h:55
@ eTestIsolatedMountain
Definition: MMFSWE.h:48
@ eTestRossbyWave
Definition: MMFSWE.h:50
const char *const TestTypeMap[]
Definition: MMFDiffusion.h:58
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
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 Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.hpp:352
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 Vvtvm(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)
vvtvm (vector times vector minus vector): z = w*x - y
Definition: Vmath.hpp:381
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
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.hpp:126
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max.
Definition: Vmath.hpp:685
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
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 > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285