Nektar++
NonlinearSWE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NonlinearSWE.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Nonlinear Shallow water equations in conservative variables
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <iomanip>
38 #include <boost/algorithm/string.hpp>
39 
42 
43 namespace Nektar
44 {
45  string NonlinearSWE::className =
47  "NonlinearSWE", NonlinearSWE::create,
48  "Nonlinear shallow water equation in conservative variables.");
49 
52  : ShallowWaterSystem(pSession)
53  {
54  }
55 
57  {
59 
61  {
64  }
65  else
66  {
67  ASSERTL0(false, "Implicit SWE not set up.");
68  }
69 
70  // Type of advection class to be used
71  switch(m_projectionType)
72  {
73  // Continuous field
75  {
76  // Do nothing
77  break;
78  }
79  // Discontinuous field
81  {
82  string advName;
83  string diffName;
84  string riemName;
85 
86  //---------------------------------------------------------------
87  // Setting up advection and diffusion operators
88  // NB: diffusion not set up for SWE at the moment
89  // but kept here for future use ...
90  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
91  // m_session->LoadSolverInfo("DiffusionType", diffName, "LDGEddy");
93  .CreateInstance(advName, advName);
94  // m_diffusion = SolverUtils::GetDiffusionFactory()
95  // .CreateInstance(diffName, diffName);
96 
97  m_advection->SetFluxVector(&NonlinearSWE::GetFluxVector, this);
98  // m_diffusion->SetFluxVectorNS(&ShallowWaterSystem::
99  // GetEddyViscosityFluxVector, this);
100 
101  // Setting up Riemann solver for advection operator
102  m_session->LoadSolverInfo("UpwindType", riemName, "Average");
104  .CreateInstance(riemName);
105 
106  // Setting up upwind solver for diffusion operator
107  // m_riemannSolverLDG = SolverUtils::GetRiemannSolverFactory()
108  // .CreateInstance("UpwindLDG");
109 
110  // Setting up parameters for advection operator Riemann solver
111  m_riemannSolver->SetParam (
112  "gravity",
113  &NonlinearSWE::GetGravity, this);
114  m_riemannSolver->SetAuxVec(
115  "vecLocs",
116  &NonlinearSWE::GetVecLocs, this);
117  m_riemannSolver->SetVector(
118  "N",
119  &NonlinearSWE::GetNormals, this);
120  m_riemannSolver->SetScalar(
121  "depth",
122  &NonlinearSWE::GetDepth, this);
123 
124  // Setting up parameters for diffusion operator Riemann solver
125  // m_riemannSolverLDG->AddParam (
126  // "gravity",
127  // &NonlinearSWE::GetGravity, this);
128  // m_riemannSolverLDG->SetAuxVec(
129  // "vecLocs",
130  // &NonlinearSWE::GetVecLocs, this);
131  // m_riemannSolverLDG->AddVector(
132  // "N",
133  // &NonlinearSWE::GetNormals, this);
134 
135  // Concluding initialisation of advection / diffusion operators
136  m_advection->SetRiemannSolver (m_riemannSolver);
137  //m_diffusion->SetRiemannSolver (m_riemannSolverLDG);
138  m_advection->InitObject (m_session, m_fields);
139  //m_diffusion->InitObject (m_session, m_fields);
140  break;
141  }
142  default:
143  {
144  ASSERTL0(false, "Unsupported projection type.");
145  break;
146  }
147  }
148 
149 
150  }
151 
153  {
154 
155  }
156 
157  // physarray contains the conservative variables
159  Array<OneD, Array<OneD, NekDouble> > &outarray)
160  {
161 
162  int ncoeffs = GetNcoeffs();
163  int nq = GetTotPoints();
164 
165  Array<OneD, NekDouble> tmp(nq);
166  Array<OneD, NekDouble> mod(ncoeffs);
167 
168  switch(m_projectionType)
169  {
171  {
172  // add to hu equation
173  Vmath::Vmul(nq,m_coriolis,1,physarray[2],1,tmp,1);
174  m_fields[0]->IProductWRTBase(tmp,mod);
175  m_fields[0]->MultiplyByElmtInvMass(mod,mod);
176  m_fields[0]->BwdTrans(mod,tmp);
177  Vmath::Vadd(nq,tmp,1,outarray[1],1,outarray[1],1);
178 
179  // add to hv equation
180  Vmath::Vmul(nq,m_coriolis,1,physarray[1],1,tmp,1);
181  Vmath::Neg(nq,tmp,1);
182  m_fields[0]->IProductWRTBase(tmp,mod);
183  m_fields[0]->MultiplyByElmtInvMass(mod,mod);
184  m_fields[0]->BwdTrans(mod,tmp);
185  Vmath::Vadd(nq,tmp,1,outarray[2],1,outarray[2],1);
186  }
187  break;
190  {
191  // add to hu equation
192  Vmath::Vmul(nq,m_coriolis,1,physarray[2],1,tmp,1);
193  Vmath::Vadd(nq,tmp,1,outarray[1],1,outarray[1],1);
194 
195  // add to hv equation
196  Vmath::Vmul(nq,m_coriolis,1,physarray[1],1,tmp,1);
197  Vmath::Neg(nq,tmp,1);
198  Vmath::Vadd(nq,tmp,1,outarray[2],1,outarray[2],1);
199  }
200  break;
201  default:
202  ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
203  break;
204  }
205 
206 
207  }
208 
209 
210  // physarray contains the conservative variables
212  Array<OneD, Array<OneD, NekDouble> > &outarray)
213  {
214 
215  int ncoeffs = GetNcoeffs();
216  int nq = GetTotPoints();
217 
218  Array<OneD, NekDouble> tmp(nq);
219  Array<OneD, NekDouble> mod(ncoeffs);
220 
221  switch(m_projectionType)
222  {
224  {
225  for (int i = 0; i < m_spacedim; ++i)
226  {
227  Vmath::Vmul(nq,m_bottomSlope[i],1,physarray[0],1,tmp,1);
228  Vmath::Smul(nq,m_g,tmp,1,tmp,1);
229  m_fields[0]->IProductWRTBase(tmp,mod);
230  m_fields[0]->MultiplyByElmtInvMass(mod,mod);
231  m_fields[0]->BwdTrans(mod,tmp);
232  Vmath::Vadd(nq,tmp,1,outarray[i+1],1,outarray[i+1],1);
233  }
234  }
235  break;
238  {
239  for (int i = 0; i < m_spacedim; ++i)
240  {
241  Vmath::Vmul(nq,m_bottomSlope[i],1,physarray[0],1,tmp,1);
242  Vmath::Smul(nq,m_g,tmp,1,tmp,1);
243  Vmath::Vadd(nq,tmp,1,outarray[i+1],1,outarray[i+1],1);
244  }
245  }
246  break;
247  default:
248  ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
249  break;
250  }
251 
252 
253  }
254 
256  Array<OneD, Array<OneD, NekDouble> >&outarray,
257  const NekDouble time)
258  {
259  int i, j;
260  int ndim = m_spacedim;
261  int nvariables = inarray.num_elements();
262  int nq = GetTotPoints();
263 
264 
265  switch(m_projectionType)
266  {
268  {
269 
270  //-------------------------------------------------------
271  // Compute the DG advection including the numerical flux
272  // by using SolverUtils/Advection
273  // Input and output in physical space
275 
276  m_advection->Advect(nvariables, m_fields, advVel, inarray,
277  outarray, time);
278  //-------------------------------------------------------
279 
280 
281  //-------------------------------------------------------
282  // negate the outarray since moving terms to the rhs
283  for(i = 0; i < nvariables; ++i)
284  {
285  Vmath::Neg(nq,outarray[i],1);
286  }
287  //-------------------------------------------------------
288 
289 
290  //-------------------------------------------------
291  // Add "source terms"
292  // Input and output in physical space
293 
294  // Coriolis forcing
295  if (m_coriolis.num_elements() != 0)
296  {
297  AddCoriolis(inarray,outarray);
298  }
299 
300  // Variable Depth
301  if (m_constantDepth != true)
302  {
303  AddVariableDepth(inarray,outarray);
304  }
305  //-------------------------------------------------
306 
307  }
308  break;
311  {
312 
313  //-------------------------------------------------------
314  // Compute the fluxvector in physical space
316  fluxvector(nvariables);
317 
318  for (i = 0; i < nvariables; ++i)
319  {
320  fluxvector[i] = Array<OneD, Array<OneD, NekDouble> >(ndim);
321  for(j = 0; j < ndim; ++j)
322  {
323  fluxvector[i][j] = Array<OneD, NekDouble>(nq);
324  }
325  }
326 
327  NonlinearSWE::GetFluxVector(inarray, fluxvector);
328  //-------------------------------------------------------
329 
330 
331 
332  //-------------------------------------------------------
333  // Take the derivative of the flux terms
334  // and negate the outarray since moving terms to the rhs
335  Array<OneD,NekDouble> tmp(nq);
336  Array<OneD, NekDouble>tmp1(nq);
337 
338  for(i = 0; i < nvariables; ++i)
339  {
340  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[0],fluxvector[i][0],tmp);
341  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[1],fluxvector[i][1],tmp1);
342  Vmath::Vadd(nq,tmp,1,tmp1,1,outarray[i],1);
343  Vmath::Neg(nq,outarray[i],1);
344  }
345 
346 
347  //-------------------------------------------------
348  // Add "source terms"
349  // Input and output in physical space
350 
351  // Coriolis forcing
352  if (m_coriolis.num_elements() != 0)
353  {
354  AddCoriolis(inarray,outarray);
355  }
356 
357  // Variable Depth
358  if (m_constantDepth != true)
359  {
360  AddVariableDepth(inarray,outarray);
361  }
362  //-------------------------------------------------
363  }
364  break;
365  default:
366  ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
367  break;
368  }
369  }
370 
371 
373  Array<OneD, Array<OneD, NekDouble> >&outarray,
374  const NekDouble time)
375  {
376  int i;
377  int nvariables = inarray.num_elements();
378 
379 
380  switch(m_projectionType)
381  {
383  {
384 
385  // Just copy over array
386  int npoints = GetNpoints();
387 
388  for(i = 0; i < nvariables; ++i)
389  {
390  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
391  }
392  SetBoundaryConditions(outarray, time);
393  break;
394  }
397  {
398 
399  EquationSystem::SetBoundaryConditions(time);
401 
402  for(i = 0; i < nvariables; ++i)
403  {
404  m_fields[i]->FwdTrans(inarray[i],coeffs);
405  m_fields[i]->BwdTrans_IterPerExp(coeffs,outarray[i]);
406  }
407  break;
408  }
409  default:
410  ASSERTL0(false,"Unknown projection scheme");
411  break;
412  }
413  }
414 
415 
416  //----------------------------------------------------
418  Array<OneD, Array<OneD, NekDouble> > &inarray,
419  NekDouble time)
420  {
421  std::string varName;
422  int nvariables = m_fields.num_elements();
423  int cnt = 0;
424 
425  // Loop over Boundary Regions
426  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
427  {
428 
429  // Wall Boundary Condition
430  if (boost::iequals(m_fields[0]->GetBndConditions()[n]->GetUserDefined(),"Wall"))
431  {
432  WallBoundary2D(n, cnt, inarray);
433  }
434 
435  // Time Dependent Boundary Condition (specified in meshfile)
436  if (m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
437  {
438  for (int i = 0; i < nvariables; ++i)
439  {
440  varName = m_session->GetVariable(i);
441  m_fields[i]->EvaluateBoundaryConditions(time, varName);
442  }
443  }
444  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
445  }
446  }
447 
448  //----------------------------------------------------
449  /**
450  * @brief Wall boundary condition.
451  */
453  int bcRegion,
454  int cnt,
455  Array<OneD, Array<OneD, NekDouble> > &physarray)
456  {
457  int i;
458  int nTracePts = GetTraceTotPoints();
459  int nvariables = physarray.num_elements();
460 
461  // get physical values of the forward trace
462  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
463  for (i = 0; i < nvariables; ++i)
464  {
465  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
466  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
467  }
468 
469  // Adjust the physical values of the trace to take
470  // user defined boundaries into account
471  int e, id1, id2, npts;
472 
473  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]
474  ->GetExpSize(); ++e)
475  {
476  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
477  GetExp(e)->GetTotPoints();
478  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
479  GetPhys_Offset(e);
480  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
481  m_fields[0]->GetTraceMap()->
482  GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
483 
484  // For 2D/3D, define: v* = v - 2(v.n)n
485  Array<OneD, NekDouble> tmp(npts, 0.0);
486 
487  // Calculate (v.n)
488  for (i = 0; i < m_spacedim; ++i)
489  {
490  Vmath::Vvtvp(npts,
491  &Fwd[1+i][id2], 1,
492  &m_traceNormals[i][id2], 1,
493  &tmp[0], 1,
494  &tmp[0], 1);
495  }
496 
497  // Calculate 2.0(v.n)
498  Vmath::Smul(npts, -2.0, &tmp[0], 1, &tmp[0], 1);
499 
500  // Calculate v* = v - 2.0(v.n)n
501  for (i = 0; i < m_spacedim; ++i)
502  {
503  Vmath::Vvtvp(npts,
504  &tmp[0], 1,
505  &m_traceNormals[i][id2], 1,
506  &Fwd[1+i][id2], 1,
507  &Fwd[1+i][id2], 1);
508  }
509 
510  // copy boundary adjusted values into the boundary expansion
511  for (i = 0; i < nvariables; ++i)
512  {
513  Vmath::Vcopy(npts, &Fwd[i][id2], 1,
514  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
515  UpdatePhys())[id1], 1);
516  }
517  }
518  }
519 
520 
521  void NonlinearSWE::WallBoundary2D(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &physarray)
522  {
523 
524  int i;
525  int nTraceNumPoints = GetTraceTotPoints();
526  int nvariables = physarray.num_elements();
527 
528  // get physical values of the forward trace
529  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
530  for (i = 0; i < nvariables; ++i)
531  {
532  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
533  m_fields[i]->ExtractTracePhys(physarray[i],Fwd[i]);
534  }
535 
536  // Adjust the physical values of the trace to take
537  // user defined boundaries into account
538  int e, id1, id2, npts;
539 
540  for(e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize(); ++e)
541  {
542  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExp(e)->GetNumPoints(0);
543  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e) ;
544  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
545 
546  switch(m_expdim)
547  {
548  case 1:
549  {
550  // negate the forward flux
551  Vmath::Neg(npts,&Fwd[1][id2],1);
552  }
553  break;
554  case 2:
555  {
556  Array<OneD, NekDouble> tmp_n(npts);
557  Array<OneD, NekDouble> tmp_t(npts);
558 
559  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[0][id2],1,&tmp_n[0],1);
560  Vmath::Vvtvp(npts,&Fwd[2][id2],1,&m_traceNormals[1][id2],1,&tmp_n[0],1,&tmp_n[0],1);
561 
562  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[1][id2],1,&tmp_t[0],1);
563  Vmath::Vvtvm(npts,&Fwd[2][id2],1,&m_traceNormals[0][id2],1,&tmp_t[0],1,&tmp_t[0],1);
564 
565  // negate the normal flux
566  Vmath::Neg(npts,tmp_n,1);
567 
568  // rotate back to Cartesian
569  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[1][id2],1,&Fwd[1][id2],1);
570  Vmath::Vvtvm(npts,&tmp_n[0],1,&m_traceNormals[0][id2],1,&Fwd[1][id2],1,&Fwd[1][id2],1);
571 
572  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[0][id2],1,&Fwd[2][id2],1);
573  Vmath::Vvtvp(npts,&tmp_n[0],1,&m_traceNormals[1][id2],1,&Fwd[2][id2],1,&Fwd[2][id2],1);
574  }
575  break;
576  case 3:
577  ASSERTL0(false,"3D not implemented for Shallow Water Equations");
578  break;
579  default:
580  ASSERTL0(false,"Illegal expansion dimension");
581  }
582 
583 
584 
585  // copy boundary adjusted values into the boundary expansion
586  for (i = 0; i < nvariables; ++i)
587  {
588  Vmath::Vcopy(npts,&Fwd[i][id2], 1,&(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1],1);
589  }
590  }
591  }
592 
593 
594  // Physfield in conservative Form
596  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
598  {
599  int i, j;
600  int nq = m_fields[0]->GetTotPoints();
601 
602  NekDouble g = m_g;
604 
605  // Flux vector for the mass equation
606  for (i = 0; i < m_spacedim; ++i)
607  {
608  velocity[i] = Array<OneD, NekDouble>(nq);
609  Vmath::Vcopy(nq, physfield[i+1], 1, flux[0][i], 1);
610  }
611 
612  GetVelocityVector(physfield, velocity);
613 
614  // Put (0.5 g h h) in tmp
615  Array<OneD, NekDouble> tmp(nq);
616  Vmath::Vmul(nq, physfield[0], 1, physfield[0], 1, tmp, 1);
617  Vmath::Smul(nq, 0.5*g, tmp, 1, tmp, 1);
618 
619  // Flux vector for the momentum equations
620  for (i = 0; i < m_spacedim; ++i)
621  {
622  for (j = 0; j < m_spacedim; ++j)
623  {
624  Vmath::Vmul(nq, velocity[j], 1, physfield[i+1], 1,
625  flux[i+1][j], 1);
626  }
627 
628  // Add (0.5 g h h) to appropriate field
629  Vmath::Vadd(nq, flux[i+1][i], 1, tmp, 1, flux[i+1][i], 1);
630  }
631 
632  }
633 
635  Array<OneD, Array<OneD, NekDouble> >&physout)
636  {
637  int nq = GetTotPoints();
638 
639  if(physin.get() == physout.get())
640  {
641  // copy indata and work with tmp array
643  for (int i = 0; i < 3; ++i)
644  {
645  // deep copy
646  tmp[i] = Array<OneD, NekDouble>(nq);
647  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
648  }
649 
650  // \eta = h - d
651  Vmath::Vsub(nq,tmp[0],1,m_depth,1,physout[0],1);
652 
653  // u = hu/h
654  Vmath::Vdiv(nq,tmp[1],1,tmp[0],1,physout[1],1);
655 
656  // v = hv/ v
657  Vmath::Vdiv(nq,tmp[2],1,tmp[0],1,physout[2],1);
658  }
659  else
660  {
661  // \eta = h - d
662  Vmath::Vsub(nq,physin[0],1,m_depth,1,physout[0],1);
663 
664  // u = hu/h
665  Vmath::Vdiv(nq,physin[1],1,physin[0],1,physout[1],1);
666 
667  // v = hv/ v
668  Vmath::Vdiv(nq,physin[2],1,physin[0],1,physout[2],1);
669  }
670  }
671 
672 
674  {
675  int nq = GetTotPoints();
676 
677  // u = hu/h
678  Vmath::Vdiv(nq,m_fields[1]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
679 
680  // v = hv/ v
681  Vmath::Vdiv(nq,m_fields[2]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
682 
683  // \eta = h - d
684  Vmath::Vsub(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
685  }
686 
688  Array<OneD, Array<OneD, NekDouble> >&physout)
689  {
690 
691  int nq = GetTotPoints();
692 
693  if(physin.get() == physout.get())
694  {
695  // copy indata and work with tmp array
697  for (int i = 0; i < 3; ++i)
698  {
699  // deep copy
700  tmp[i] = Array<OneD, NekDouble>(nq);
701  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
702  }
703 
704  // h = \eta + d
705  Vmath::Vadd(nq,tmp[0],1,m_depth,1,physout[0],1);
706 
707  // hu = h * u
708  Vmath::Vmul(nq,physout[0],1,tmp[1],1,physout[1],1);
709 
710  // hv = h * v
711  Vmath::Vmul(nq,physout[0],1,tmp[2],1,physout[2],1);
712 
713  }
714  else
715  {
716  // h = \eta + d
717  Vmath::Vadd(nq,physin[0],1,m_depth,1,physout[0],1);
718 
719  // hu = h * u
720  Vmath::Vmul(nq,physout[0],1,physin[1],1,physout[1],1);
721 
722  // hv = h * v
723  Vmath::Vmul(nq,physout[0],1,physin[2],1,physout[2],1);
724 
725  }
726 
727  }
728 
730  {
731  int nq = GetTotPoints();
732 
733  // h = \eta + d
734  Vmath::Vadd(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
735 
736  // hu = h * u
737  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[1]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
738 
739  // hv = h * v
740  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[2]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
741  }
742 
743 
744  /**
745  * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
746  * \f$ h\mathbf{v} \f$.
747  *
748  * @param physfield Momentum field.
749  * @param velocity Velocity field.
750  */
752  const Array<OneD, Array<OneD, NekDouble> > &physfield,
753  Array<OneD, Array<OneD, NekDouble> > &velocity)
754  {
755  const int npts = physfield[0].num_elements();
756 
757  for (int i = 0; i < m_spacedim; ++i)
758  {
759  Vmath::Vdiv(npts, physfield[1+i], 1, physfield[0], 1,
760  velocity[i], 1);
761  }
762  }
763 
764 
766  {
768  SolverUtils::AddSummaryItem(s, "Variables", "h should be in field[0]");
769  SolverUtils::AddSummaryItem(s, "", "hu should be in field[1]");
770  SolverUtils::AddSummaryItem(s, "", "hv should be in field[2]");
771  }
772 
773 } //end of namespace
774 
Array< OneD, NekDouble > m_coriolis
Coriolis force.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
virtual void v_PrimitiveToConservative()
Base class for unsteady solvers.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Array< OneD, NekDouble > m_depth
Still water depth.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
int m_expdim
Expansion dimension.
static std::string className
Name of class.
Definition: NonlinearSWE.h:64
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:428
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
SolverUtils::AdvectionSharedPtr m_advection
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
Definition: NonlinearSWE.h:56
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:227
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
void AddVariableDepth(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
const Array< OneD, NekDouble > & GetDepth()
virtual void v_InitObject()
Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT int GetTotPoints()
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void WallBoundary(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary condition.
Array< OneD, Array< OneD, NekDouble > > m_bottomSlope
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual void v_InitObject()
Init object for UnsteadySystem class.
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
RiemannSolverFactory & GetRiemannSolverFactory()
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
static std::string npts
Definition: InputFld.cpp:43
int m_spacedim
Spatial dimension (>= expansion dim).
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
EquationSystemFactory & GetEquationSystemFactory()
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:329
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
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:451
virtual void v_ConservativeToPrimitive()
SOLVER_UTILS_EXPORT int GetNcoeffs()
NonlinearSWE(const LibUtilities::SessionReaderSharedPtr &pSession)
bool m_constantDepth
Indicates if constant depth case.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
NekDouble m_g
Acceleration of gravity.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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:285
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:169
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215