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