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