Nektar++
LinearSWE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File LinearSWE.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: Linear Shallow water equations in primitive variables
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <iomanip>
38 #include <boost/algorithm/string.hpp>
39 
42 
43 namespace Nektar
44 {
45  string LinearSWE::className =
47  "LinearSWE", LinearSWE::create,
48  "Linear shallow water equation in primitive 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(&LinearSWE::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, "NoSolver");
103  if ((riemName == "LinearHLL") && (m_constantDepth != true))
104  {
105  ASSERTL0(false,"LinearHLL only valid for constant depth");
106  }
108  .CreateInstance(riemName);
109 
110  // Setting up upwind solver for diffusion operator
111  // m_riemannSolverLDG = SolverUtils::GetRiemannSolverFactory()
112  // .CreateInstance("UpwindLDG");
113 
114  // Setting up parameters for advection operator Riemann solver
115  m_riemannSolver->SetParam (
116  "gravity",
117  &LinearSWE::GetGravity, this);
118  m_riemannSolver->SetAuxVec(
119  "vecLocs",
120  &LinearSWE::GetVecLocs, this);
121  m_riemannSolver->SetVector(
122  "N",
123  &LinearSWE::GetNormals, this);
124 
125  // The numerical flux for linear SWE requires depth information
126  int nTracePointsTot = m_fields[0]->GetTrace()->GetTotPoints();
127  m_dFwd = Array<OneD, NekDouble>(nTracePointsTot);
128  m_dBwd = Array<OneD, NekDouble>(nTracePointsTot);
129  m_fields[0]->GetFwdBwdTracePhys(m_depth, m_dFwd, m_dBwd);
130  CopyBoundaryTrace(m_dFwd,m_dBwd);
131  m_riemannSolver->SetScalar(
132  "depthFwd",
133  &LinearSWE::GetDepthFwd, this);
134  m_riemannSolver->SetScalar(
135  "depthBwd",
136  &LinearSWE::GetDepthBwd, this);
137 
138  // Setting up parameters for diffusion operator Riemann solver
139  // m_riemannSolverLDG->AddParam (
140  // "gravity",
141  // &NonlinearSWE::GetGravity, this);
142  // m_riemannSolverLDG->SetAuxVec(
143  // "vecLocs",
144  // &NonlinearSWE::GetVecLocs, this);
145  // m_riemannSolverLDG->AddVector(
146  // "N",
147  // &NonlinearSWE::GetNormals, this);
148 
149  // Concluding initialisation of advection / diffusion operators
150  m_advection->SetRiemannSolver (m_riemannSolver);
151  //m_diffusion->SetRiemannSolver (m_riemannSolverLDG);
152  m_advection->InitObject (m_session, m_fields);
153  //m_diffusion->InitObject (m_session, m_fields);
154  break;
155  }
156  default:
157  {
158  ASSERTL0(false, "Unsupported projection type.");
159  break;
160  }
161  }
162 
163 
164  }
165 
167  {
168 
169  }
170 
171  // physarray contains the conservative variables
172  void LinearSWE::AddCoriolis(const Array<OneD, const Array<OneD, NekDouble> > &physarray,
173  Array<OneD, Array<OneD, NekDouble> > &outarray)
174  {
175 
176  int ncoeffs = GetNcoeffs();
177  int nq = GetTotPoints();
178 
179  Array<OneD, NekDouble> tmp(nq);
180  Array<OneD, NekDouble> mod(ncoeffs);
181 
182  switch(m_projectionType)
183  {
185  {
186  // add to u equation
187  Vmath::Vmul(nq,m_coriolis,1,physarray[2],1,tmp,1);
188  m_fields[0]->IProductWRTBase(tmp,mod);
189  m_fields[0]->MultiplyByElmtInvMass(mod,mod);
190  m_fields[0]->BwdTrans(mod,tmp);
191  Vmath::Vadd(nq,tmp,1,outarray[1],1,outarray[1],1);
192 
193  // add to v equation
194  Vmath::Vmul(nq,m_coriolis,1,physarray[1],1,tmp,1);
195  Vmath::Neg(nq,tmp,1);
196  m_fields[0]->IProductWRTBase(tmp,mod);
197  m_fields[0]->MultiplyByElmtInvMass(mod,mod);
198  m_fields[0]->BwdTrans(mod,tmp);
199  Vmath::Vadd(nq,tmp,1,outarray[2],1,outarray[2],1);
200  }
201  break;
204  {
205  // add to u equation
206  Vmath::Vmul(nq,m_coriolis,1,physarray[2],1,tmp,1);
207  Vmath::Vadd(nq,tmp,1,outarray[1],1,outarray[1],1);
208 
209  // add to v equation
210  Vmath::Vmul(nq,m_coriolis,1,physarray[1],1,tmp,1);
211  Vmath::Neg(nq,tmp,1);
212  Vmath::Vadd(nq,tmp,1,outarray[2],1,outarray[2],1);
213  }
214  break;
215  default:
216  ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
217  break;
218  }
219 
220 
221  }
222 
223  void LinearSWE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
224  Array<OneD, Array<OneD, NekDouble> >&outarray,
225  const NekDouble time)
226  {
227  int i, j;
228  int ndim = m_spacedim;
229  int nvariables = inarray.num_elements();
230  int nq = GetTotPoints();
231 
232 
233  switch(m_projectionType)
234  {
236  {
237 
238  //-------------------------------------------------------
239  // Compute the DG advection including the numerical flux
240  // by using SolverUtils/Advection
241  // Input and output in physical space
243 
244  m_advection->Advect(nvariables, m_fields, advVel, inarray,
245  outarray, time);
246  //-------------------------------------------------------
247 
248 
249  //-------------------------------------------------------
250  // negate the outarray since moving terms to the rhs
251  for(i = 0; i < nvariables; ++i)
252  {
253  Vmath::Neg(nq,outarray[i],1);
254  }
255  //-------------------------------------------------------
256 
257 
258  //-------------------------------------------------
259  // Add "source terms"
260  // Input and output in physical space
261 
262  // Coriolis forcing
263  if (m_coriolis.num_elements() != 0)
264  {
265  AddCoriolis(inarray,outarray);
266  }
267  //-------------------------------------------------
268 
269  }
270  break;
273  {
274 
275  //-------------------------------------------------------
276  // Compute the fluxvector in physical space
278  fluxvector(nvariables);
279 
280  for (i = 0; i < nvariables; ++i)
281  {
282  fluxvector[i] = Array<OneD, Array<OneD, NekDouble> >(ndim);
283  for(j = 0; j < ndim; ++j)
284  {
285  fluxvector[i][j] = Array<OneD, NekDouble>(nq);
286  }
287  }
288 
289  LinearSWE::GetFluxVector(inarray, fluxvector);
290  //-------------------------------------------------------
291 
292 
293  //-------------------------------------------------------
294  // Take the derivative of the flux terms
295  // and negate the outarray since moving terms to the rhs
296  Array<OneD,NekDouble> tmp(nq);
297  Array<OneD, NekDouble>tmp1(nq);
298 
299  for(i = 0; i < nvariables; ++i)
300  {
301  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[0],fluxvector[i][0],tmp);
302  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[1],fluxvector[i][1],tmp1);
303  Vmath::Vadd(nq,tmp,1,tmp1,1,outarray[i],1);
304  Vmath::Neg(nq,outarray[i],1);
305  }
306 
307  //-------------------------------------------------
308  // Add "source terms"
309  // Input and output in physical space
310 
311  // Coriolis forcing
312  if (m_coriolis.num_elements() != 0)
313  {
314  AddCoriolis(inarray,outarray);
315  }
316  //-------------------------------------------------
317  }
318  break;
319  default:
320  ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
321  break;
322  }
323  }
324 
325 
327  Array<OneD, Array<OneD, NekDouble> >&outarray,
328  const NekDouble time)
329  {
330  int i;
331  int nvariables = inarray.num_elements();
332 
333 
334  switch(m_projectionType)
335  {
337  {
338 
339  // Just copy over array
340  int npoints = GetNpoints();
341 
342  for(i = 0; i < nvariables; ++i)
343  {
344  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
345  }
346  SetBoundaryConditions(outarray, time);
347  break;
348  }
351  {
352 
353  EquationSystem::SetBoundaryConditions(time);
355 
356  for(i = 0; i < nvariables; ++i)
357  {
358  m_fields[i]->FwdTrans(inarray[i],coeffs);
359  m_fields[i]->BwdTrans_IterPerExp(coeffs,outarray[i]);
360  }
361  break;
362  }
363  default:
364  ASSERTL0(false,"Unknown projection scheme");
365  break;
366  }
367  }
368 
369 
370  //----------------------------------------------------
372  Array<OneD, Array<OneD, NekDouble> > &inarray,
373  NekDouble time)
374  {
375  std::string varName;
376  int nvariables = m_fields.num_elements();
377  int cnt = 0;
378 
379  // loop over Boundary Regions
380  for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
381  {
382  // Wall Boundary Condition
383  if (boost::iequals(m_fields[0]->GetBndConditions()[n]->GetUserDefined(),"Wall"))
384  {
385  WallBoundary2D(n, cnt, inarray);
386  }
387 
388  // Time Dependent Boundary Condition (specified in meshfile)
389  if (m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
390  {
391  for (int i = 0; i < nvariables; ++i)
392  {
393  varName = m_session->GetVariable(i);
394  m_fields[i]->EvaluateBoundaryConditions(time, varName);
395  }
396  }
397  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
398  }
399  }
400 
401  //----------------------------------------------------
402  /**
403  * @brief Wall boundary condition.
404  */
406  int bcRegion,
407  int cnt,
408  Array<OneD, Array<OneD, NekDouble> > &physarray)
409  {
410  int i;
411  int nTracePts = GetTraceTotPoints();
412  int nvariables = physarray.num_elements();
413 
414  // get physical values of the forward trace
415  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
416  for (i = 0; i < nvariables; ++i)
417  {
418  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
419  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
420  }
421 
422  // Adjust the physical values of the trace to take
423  // user defined boundaries into account
424  int e, id1, id2, npts;
425 
426  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]
427  ->GetExpSize(); ++e)
428  {
429  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
430  GetExp(e)->GetTotPoints();
431  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
432  GetPhys_Offset(e);
433  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
434  m_fields[0]->GetTraceMap()->
435  GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
436 
437  // For 2D/3D, define: v* = v - 2(v.n)n
438  Array<OneD, NekDouble> tmp(npts, 0.0);
439 
440  // Calculate (v.n)
441  for (i = 0; i < m_spacedim; ++i)
442  {
443  Vmath::Vvtvp(npts,
444  &Fwd[1+i][id2], 1,
445  &m_traceNormals[i][id2], 1,
446  &tmp[0], 1,
447  &tmp[0], 1);
448  }
449 
450  // Calculate 2.0(v.n)
451  Vmath::Smul(npts, -2.0, &tmp[0], 1, &tmp[0], 1);
452 
453  // Calculate v* = v - 2.0(v.n)n
454  for (i = 0; i < m_spacedim; ++i)
455  {
456  Vmath::Vvtvp(npts,
457  &tmp[0], 1,
458  &m_traceNormals[i][id2], 1,
459  &Fwd[1+i][id2], 1,
460  &Fwd[1+i][id2], 1);
461  }
462 
463  // copy boundary adjusted values into the boundary expansion
464  for (i = 0; i < nvariables; ++i)
465  {
466  Vmath::Vcopy(npts, &Fwd[i][id2], 1,
467  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
468  UpdatePhys())[id1], 1);
469  }
470  }
471  }
472 
473 
474  void LinearSWE::WallBoundary2D(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &physarray)
475  {
476 
477  int i;
478  int nTraceNumPoints = GetTraceTotPoints();
479  int nvariables = physarray.num_elements();
480 
481  // get physical values of the forward trace
482  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
483  for (i = 0; i < nvariables; ++i)
484  {
485  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
486  m_fields[i]->ExtractTracePhys(physarray[i],Fwd[i]);
487  }
488 
489  // Adjust the physical values of the trace to take
490  // user defined boundaries into account
491  int e, id1, id2, npts;
492 
493  for(e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize(); ++e)
494  {
495  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExp(e)->GetNumPoints(0);
496  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e) ;
497  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
498 
499  switch(m_expdim)
500  {
501  case 1:
502  {
503  // negate the forward flux
504  Vmath::Neg(npts,&Fwd[1][id2],1);
505  }
506  break;
507  case 2:
508  {
509  Array<OneD, NekDouble> tmp_n(npts);
510  Array<OneD, NekDouble> tmp_t(npts);
511 
512  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[0][id2],1,&tmp_n[0],1);
513  Vmath::Vvtvp(npts,&Fwd[2][id2],1,&m_traceNormals[1][id2],1,&tmp_n[0],1,&tmp_n[0],1);
514 
515  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[1][id2],1,&tmp_t[0],1);
516  Vmath::Vvtvm(npts,&Fwd[2][id2],1,&m_traceNormals[0][id2],1,&tmp_t[0],1,&tmp_t[0],1);
517 
518  // negate the normal flux
519  Vmath::Neg(npts,tmp_n,1);
520 
521  // rotate back to Cartesian
522  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[1][id2],1,&Fwd[1][id2],1);
523  Vmath::Vvtvm(npts,&tmp_n[0],1,&m_traceNormals[0][id2],1,&Fwd[1][id2],1,&Fwd[1][id2],1);
524 
525  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[0][id2],1,&Fwd[2][id2],1);
526  Vmath::Vvtvp(npts,&tmp_n[0],1,&m_traceNormals[1][id2],1,&Fwd[2][id2],1,&Fwd[2][id2],1);
527  }
528  break;
529  case 3:
530  ASSERTL0(false,"3D not implemented for Shallow Water Equations");
531  break;
532  default:
533  ASSERTL0(false,"Illegal expansion dimension");
534  }
535 
536 
537 
538  // copy boundary adjusted values into the boundary expansion
539  for (i = 0; i < nvariables; ++i)
540  {
541  Vmath::Vcopy(npts,&Fwd[i][id2], 1,&(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1],1);
542  }
543  }
544  }
545 
546 
547  // Physfield in primitive Form
549  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
551  {
552  int i, j;
553  int nq = m_fields[0]->GetTotPoints();
554 
555  NekDouble g = m_g;
556 
557  // Flux vector for the mass equation
558  for (i = 0; i < m_spacedim; ++i)
559  {
560  Vmath::Vmul(nq, m_depth, 1, physfield[i+1], 1, flux[0][i], 1);
561  }
562 
563  // Put (g eta) in tmp
564  Array<OneD, NekDouble> tmp(nq);
565  Vmath::Smul(nq, g, physfield[0], 1, tmp, 1);
566 
567  // Flux vector for the momentum equations
568  for (i = 0; i < m_spacedim; ++i)
569  {
570  for (j = 0; j < m_spacedim; ++j)
571  {
572  // must zero fluxes as not initialised to zero in AdvectionWeakDG ...
573  Vmath::Zero(nq, flux[i+1][j], 1);
574  }
575 
576  // Add (g eta) to appropriate field
577  Vmath::Vadd(nq, flux[i+1][i], 1, tmp, 1, flux[i+1][i], 1);
578  }
579 
580  }
581 
583  Array<OneD, Array<OneD, NekDouble> >&physout)
584  {
585  int nq = GetTotPoints();
586 
587  if(physin.get() == physout.get())
588  {
589  // copy indata and work with tmp array
591  for (int i = 0; i < 3; ++i)
592  {
593  // deep copy
594  tmp[i] = Array<OneD, NekDouble>(nq);
595  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
596  }
597 
598  // \eta = h - d
599  Vmath::Vsub(nq,tmp[0],1,m_depth,1,physout[0],1);
600 
601  // u = hu/h
602  Vmath::Vdiv(nq,tmp[1],1,tmp[0],1,physout[1],1);
603 
604  // v = hv/ v
605  Vmath::Vdiv(nq,tmp[2],1,tmp[0],1,physout[2],1);
606  }
607  else
608  {
609  // \eta = h - d
610  Vmath::Vsub(nq,physin[0],1,m_depth,1,physout[0],1);
611 
612  // u = hu/h
613  Vmath::Vdiv(nq,physin[1],1,physin[0],1,physout[1],1);
614 
615  // v = hv/ v
616  Vmath::Vdiv(nq,physin[2],1,physin[0],1,physout[2],1);
617  }
618  }
619 
620 
622  {
623  int nq = GetTotPoints();
624 
625  // u = hu/h
626  Vmath::Vdiv(nq,m_fields[1]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
627 
628  // v = hv/ v
629  Vmath::Vdiv(nq,m_fields[2]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
630 
631  // \eta = h - d
632  Vmath::Vsub(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
633  }
634 
636  Array<OneD, Array<OneD, NekDouble> >&physout)
637  {
638 
639  int nq = GetTotPoints();
640 
641  if(physin.get() == physout.get())
642  {
643  // copy indata and work with tmp array
645  for (int i = 0; i < 3; ++i)
646  {
647  // deep copy
648  tmp[i] = Array<OneD, NekDouble>(nq);
649  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
650  }
651 
652  // h = \eta + d
653  Vmath::Vadd(nq,tmp[0],1,m_depth,1,physout[0],1);
654 
655  // hu = h * u
656  Vmath::Vmul(nq,physout[0],1,tmp[1],1,physout[1],1);
657 
658  // hv = h * v
659  Vmath::Vmul(nq,physout[0],1,tmp[2],1,physout[2],1);
660 
661  }
662  else
663  {
664  // h = \eta + d
665  Vmath::Vadd(nq,physin[0],1,m_depth,1,physout[0],1);
666 
667  // hu = h * u
668  Vmath::Vmul(nq,physout[0],1,physin[1],1,physout[1],1);
669 
670  // hv = h * v
671  Vmath::Vmul(nq,physout[0],1,physin[2],1,physout[2],1);
672 
673  }
674 
675  }
676 
678  {
679  int nq = GetTotPoints();
680 
681  // h = \eta + d
682  Vmath::Vadd(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
683 
684  // hu = h * u
685  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[1]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
686 
687  // hv = h * v
688  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[2]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
689  }
690 
691 
692  /**
693  * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
694  * \f$ h\mathbf{v} \f$.
695  *
696  * @param physfield Velocity field.
697  * @param velocity Velocity field.
698  */
700  const Array<OneD, Array<OneD, NekDouble> > &physfield,
701  Array<OneD, Array<OneD, NekDouble> > &velocity)
702  {
703  const int npts = physfield[0].num_elements();
704 
705  for (int i = 0; i < m_spacedim; ++i)
706  {
707  Vmath::Vcopy(npts, physfield[1+i], 1, velocity[i], 1);
708  }
709  }
710 
711 
713  {
715  if (m_session->DefinesSolverInfo("UpwindType"))
716  {
717  std::string UpwindType;
718  UpwindType = m_session->GetSolverInfo("UpwindType");
719  if (UpwindType == "LinearAverage")
720  {
721  SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear Average");
722  }
723  if (UpwindType == "LinearHLL")
724  {
725  SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear HLL");
726  }
727  }
728  SolverUtils::AddSummaryItem(s, "Variables", "eta should be in field[0]");
729  SolverUtils::AddSummaryItem(s, "", "u should be in field[1]");
730  SolverUtils::AddSummaryItem(s, "", "v should be in field[2]");
731  }
732 
733 } //end of namespace
734 
Array< OneD, NekDouble > m_coriolis
Coriolis force.
const Array< OneD, NekDouble > & GetDepthFwd()
Definition: LinearSWE.h:96
void WallBoundary(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary condition.
Definition: LinearSWE.cpp:405
virtual void v_InitObject()
Init object for UnsteadySystem class.
Definition: LinearSWE.cpp:56
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
Definition: LinearSWE.cpp:712
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.
Array< OneD, NekDouble > m_dFwd
Still water depth traces.
Definition: LinearSWE.h:75
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Definition: LinearSWE.cpp:548
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Definition: LinearSWE.cpp:474
virtual void v_PrimitiveToConservative()
Definition: LinearSWE.cpp:677
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
int m_expdim
Expansion dimension.
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()
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: LinearSWE.cpp:172
SolverUtils::AdvectionSharedPtr m_advection
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
Definition: LinearSWE.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
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
SOLVER_UTILS_EXPORT int GetTotPoints()
virtual void v_ConservativeToPrimitive()
Definition: LinearSWE.cpp:621
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Definition: LinearSWE.cpp:371
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 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 ~LinearSWE()
Definition: LinearSWE.cpp:166
virtual void v_InitObject()
Init object for UnsteadySystem class.
static std::string className
Name of class.
Definition: LinearSWE.h:64
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).
void CopyBoundaryTrace(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
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
Array< OneD, NekDouble > m_dBwd
Definition: LinearSWE.h:76
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition: LinearSWE.cpp:326
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()
LinearSWE(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: LinearSWE.cpp:50
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
Definition: LinearSWE.cpp:699
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
SOLVER_UTILS_EXPORT int GetNcoeffs()
bool m_constantDepth
Indicates if constant depth case.
const Array< OneD, NekDouble > & GetDepthBwd()
Definition: LinearSWE.h:100
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition: LinearSWE.cpp:223
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215