Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
40 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
41 #include <ShallowWaterSolver/EquationSystems/LinearSWE.h>
42 
43 namespace Nektar
44 {
45  string LinearSWE::className =
46  SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
47  "LinearSWE", LinearSWE::create,
48  "Linear shallow water equation in primitive variables.");
49 
50  LinearSWE::LinearSWE(
51  const LibUtilities::SessionReaderSharedPtr& pSession)
52  : ShallowWaterSystem(pSession)
53  {
54  }
55 
56  void LinearSWE::v_InitObject()
57  {
58  ShallowWaterSystem::v_InitObject();
59 
60  if (m_explicitAdvection)
61  {
62  m_ode.DefineOdeRhs (&LinearSWE::DoOdeRhs, this);
63  m_ode.DefineProjection (&LinearSWE::DoOdeProjection, this);
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
74  case MultiRegions::eGalerkin:
75  {
76  // Do nothing
77  break;
78  }
79  // Discontinuous field
80  case MultiRegions::eDiscontinuous:
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");
92  m_advection = SolverUtils::GetAdvectionFactory()
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  }
107  m_riemannSolver = SolverUtils::GetRiemannSolverFactory()
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 
166  LinearSWE::~LinearSWE()
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  {
184  case MultiRegions::eDiscontinuous:
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;
202  case MultiRegions::eGalerkin:
203  case MultiRegions::eMixed_CG_Discontinuous:
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  {
235  case MultiRegions::eDiscontinuous:
236  {
237 
238  //-------------------------------------------------------
239  // Compute the DG advection including the numerical flux
240  // by using SolverUtils/Advection
241  // Input and output in physical space
242  Array<OneD, Array<OneD, NekDouble> > advVel;
243 
244  m_advection->Advect(nvariables, m_fields, advVel, inarray, outarray);
245  //-------------------------------------------------------
246 
247 
248  //-------------------------------------------------------
249  // negate the outarray since moving terms to the rhs
250  for(i = 0; i < nvariables; ++i)
251  {
252  Vmath::Neg(nq,outarray[i],1);
253  }
254  //-------------------------------------------------------
255 
256 
257  //-------------------------------------------------
258  // Add "source terms"
259  // Input and output in physical space
260 
261  // Coriolis forcing
262  if (m_coriolis.num_elements() != 0)
263  {
264  AddCoriolis(inarray,outarray);
265  }
266  //-------------------------------------------------
267 
268  }
269  break;
270  case MultiRegions::eGalerkin:
271  case MultiRegions::eMixed_CG_Discontinuous:
272  {
273 
274  //-------------------------------------------------------
275  // Compute the fluxvector in physical space
276  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
277  fluxvector(nvariables);
278 
279  for (i = 0; i < nvariables; ++i)
280  {
281  fluxvector[i] = Array<OneD, Array<OneD, NekDouble> >(ndim);
282  for(j = 0; j < ndim; ++j)
283  {
284  fluxvector[i][j] = Array<OneD, NekDouble>(nq);
285  }
286  }
287 
288  LinearSWE::GetFluxVector(inarray, fluxvector);
289  //-------------------------------------------------------
290 
291 
292  //-------------------------------------------------------
293  // Take the derivative of the flux terms
294  // and negate the outarray since moving terms to the rhs
295  Array<OneD,NekDouble> tmp(nq);
296  Array<OneD, NekDouble>tmp1(nq);
297 
298  for(i = 0; i < nvariables; ++i)
299  {
300  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[0],fluxvector[i][0],tmp);
301  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[1],fluxvector[i][1],tmp1);
302  Vmath::Vadd(nq,tmp,1,tmp1,1,outarray[i],1);
303  Vmath::Neg(nq,outarray[i],1);
304  }
305 
306  //-------------------------------------------------
307  // Add "source terms"
308  // Input and output in physical space
309 
310  // Coriolis forcing
311  if (m_coriolis.num_elements() != 0)
312  {
313  AddCoriolis(inarray,outarray);
314  }
315  //-------------------------------------------------
316  }
317  break;
318  default:
319  ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
320  break;
321  }
322  }
323 
324 
325  void LinearSWE::DoOdeProjection(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
326  Array<OneD, Array<OneD, NekDouble> >&outarray,
327  const NekDouble time)
328  {
329  int i;
330  int nvariables = inarray.num_elements();
331 
332 
333  switch(m_projectionType)
334  {
335  case MultiRegions::eDiscontinuous:
336  {
337 
338  // Just copy over array
339  int npoints = GetNpoints();
340 
341  for(i = 0; i < nvariables; ++i)
342  {
343  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
344  }
345  SetBoundaryConditions(outarray, time);
346  break;
347  }
348  case MultiRegions::eGalerkin:
349  case MultiRegions::eMixed_CG_Discontinuous:
350  {
351 
352  EquationSystem::SetBoundaryConditions(time);
353  Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
354 
355  for(i = 0; i < nvariables; ++i)
356  {
357  m_fields[i]->FwdTrans(inarray[i],coeffs);
358  m_fields[i]->BwdTrans_IterPerExp(coeffs,outarray[i]);
359  }
360  break;
361  }
362  default:
363  ASSERTL0(false,"Unknown projection scheme");
364  break;
365  }
366  }
367 
368 
369  //----------------------------------------------------
370  void LinearSWE::SetBoundaryConditions(
371  Array<OneD, Array<OneD, NekDouble> > &inarray,
372  NekDouble time)
373  {
374  std::string varName;
375  int nvariables = m_fields.num_elements();
376  int cnt = 0;
377 
378  // loop over Boundary Regions
379  for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
380  {
381  // Wall Boundary Condition
382  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
383  SpatialDomains::eWall)
384  {
385  WallBoundary2D(n, cnt, inarray);
386  }
387 
388  // Time Dependent Boundary Condition (specified in meshfile)
389  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
390  SpatialDomains::eTimeDependent)
391  {
392  for (int i = 0; i < nvariables; ++i)
393  {
394  varName = m_session->GetVariable(i);
395  m_fields[i]->EvaluateBoundaryConditions(time, varName);
396  }
397  }
398  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
399  }
400  }
401 
402  //----------------------------------------------------
403  /**
404  * @brief Wall boundary condition.
405  */
406  void LinearSWE::WallBoundary(
407  int bcRegion,
408  int cnt,
409  Array<OneD, Array<OneD, NekDouble> > &physarray)
410  {
411  int i;
412  int nTracePts = GetTraceTotPoints();
413  int nvariables = physarray.num_elements();
414 
415  // get physical values of the forward trace
416  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
417  for (i = 0; i < nvariables; ++i)
418  {
419  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
420  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
421  }
422 
423  // Adjust the physical values of the trace to take
424  // user defined boundaries into account
425  int e, id1, id2, npts;
426 
427  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]
428  ->GetExpSize(); ++e)
429  {
430  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
431  GetExp(e)->GetTotPoints();
432  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
433  GetPhys_Offset(e);
434  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
435  m_fields[0]->GetTraceMap()->
436  GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
437 
438  // For 2D/3D, define: v* = v - 2(v.n)n
439  Array<OneD, NekDouble> tmp(npts, 0.0);
440 
441  // Calculate (v.n)
442  for (i = 0; i < m_spacedim; ++i)
443  {
444  Vmath::Vvtvp(npts,
445  &Fwd[1+i][id2], 1,
446  &m_traceNormals[i][id2], 1,
447  &tmp[0], 1,
448  &tmp[0], 1);
449  }
450 
451  // Calculate 2.0(v.n)
452  Vmath::Smul(npts, -2.0, &tmp[0], 1, &tmp[0], 1);
453 
454  // Calculate v* = v - 2.0(v.n)n
455  for (i = 0; i < m_spacedim; ++i)
456  {
457  Vmath::Vvtvp(npts,
458  &tmp[0], 1,
459  &m_traceNormals[i][id2], 1,
460  &Fwd[1+i][id2], 1,
461  &Fwd[1+i][id2], 1);
462  }
463 
464  // copy boundary adjusted values into the boundary expansion
465  for (i = 0; i < nvariables; ++i)
466  {
467  Vmath::Vcopy(npts, &Fwd[i][id2], 1,
468  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
469  UpdatePhys())[id1], 1);
470  }
471  }
472  }
473 
474 
475  void LinearSWE::WallBoundary2D(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &physarray)
476  {
477 
478  int i;
479  int nTraceNumPoints = GetTraceTotPoints();
480  int nvariables = physarray.num_elements();
481 
482  // get physical values of the forward trace
483  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
484  for (i = 0; i < nvariables; ++i)
485  {
486  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
487  m_fields[i]->ExtractTracePhys(physarray[i],Fwd[i]);
488  }
489 
490  // Adjust the physical values of the trace to take
491  // user defined boundaries into account
492  int e, id1, id2, npts;
493 
494  for(e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize(); ++e)
495  {
496  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExp(e)->GetNumPoints(0);
497  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e) ;
498  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
499 
500  switch(m_expdim)
501  {
502  case 1:
503  {
504  // negate the forward flux
505  Vmath::Neg(npts,&Fwd[1][id2],1);
506  }
507  break;
508  case 2:
509  {
510  Array<OneD, NekDouble> tmp_n(npts);
511  Array<OneD, NekDouble> tmp_t(npts);
512 
513  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[0][id2],1,&tmp_n[0],1);
514  Vmath::Vvtvp(npts,&Fwd[2][id2],1,&m_traceNormals[1][id2],1,&tmp_n[0],1,&tmp_n[0],1);
515 
516  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[1][id2],1,&tmp_t[0],1);
517  Vmath::Vvtvm(npts,&Fwd[2][id2],1,&m_traceNormals[0][id2],1,&tmp_t[0],1,&tmp_t[0],1);
518 
519  // negate the normal flux
520  Vmath::Neg(npts,tmp_n,1);
521 
522  // rotate back to Cartesian
523  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[1][id2],1,&Fwd[1][id2],1);
524  Vmath::Vvtvm(npts,&tmp_n[0],1,&m_traceNormals[0][id2],1,&Fwd[1][id2],1,&Fwd[1][id2],1);
525 
526  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[0][id2],1,&Fwd[2][id2],1);
527  Vmath::Vvtvp(npts,&tmp_n[0],1,&m_traceNormals[1][id2],1,&Fwd[2][id2],1,&Fwd[2][id2],1);
528  }
529  break;
530  case 3:
531  ASSERTL0(false,"3D not implemented for Shallow Water Equations");
532  break;
533  default:
534  ASSERTL0(false,"Illegal expansion dimension");
535  }
536 
537 
538 
539  // copy boundary adjusted values into the boundary expansion
540  for (i = 0; i < nvariables; ++i)
541  {
542  Vmath::Vcopy(npts,&Fwd[i][id2], 1,&(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1],1);
543  }
544  }
545  }
546 
547 
548  // Physfield in primitive Form
549  void LinearSWE::GetFluxVector(
550  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
551  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
552  {
553  int i, j;
554  int nq = m_fields[0]->GetTotPoints();
555 
556  NekDouble g = m_g;
557 
558  // Flux vector for the mass equation
559  for (i = 0; i < m_spacedim; ++i)
560  {
561  Vmath::Vmul(nq, m_depth, 1, physfield[i+1], 1, flux[0][i], 1);
562  }
563 
564  // Put (g eta) in tmp
565  Array<OneD, NekDouble> tmp(nq);
566  Vmath::Smul(nq, g, physfield[0], 1, tmp, 1);
567 
568  // Flux vector for the momentum equations
569  for (i = 0; i < m_spacedim; ++i)
570  {
571  for (j = 0; j < m_spacedim; ++j)
572  {
573  // must zero fluxes as not initialised to zero in AdvectionWeakDG ...
574  Vmath::Zero(nq, flux[i+1][j], 1);
575  }
576 
577  // Add (g eta) to appropriate field
578  Vmath::Vadd(nq, flux[i+1][i], 1, tmp, 1, flux[i+1][i], 1);
579  }
580 
581  }
582 
583  void LinearSWE::ConservativeToPrimitive(const Array<OneD, const Array<OneD, NekDouble> >&physin,
584  Array<OneD, Array<OneD, NekDouble> >&physout)
585  {
586  int nq = GetTotPoints();
587 
588  if(physin.get() == physout.get())
589  {
590  // copy indata and work with tmp array
591  Array<OneD, Array<OneD, NekDouble> >tmp(3);
592  for (int i = 0; i < 3; ++i)
593  {
594  // deep copy
595  tmp[i] = Array<OneD, NekDouble>(nq);
596  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
597  }
598 
599  // \eta = h - d
600  Vmath::Vsub(nq,tmp[0],1,m_depth,1,physout[0],1);
601 
602  // u = hu/h
603  Vmath::Vdiv(nq,tmp[1],1,tmp[0],1,physout[1],1);
604 
605  // v = hv/ v
606  Vmath::Vdiv(nq,tmp[2],1,tmp[0],1,physout[2],1);
607  }
608  else
609  {
610  // \eta = h - d
611  Vmath::Vsub(nq,physin[0],1,m_depth,1,physout[0],1);
612 
613  // u = hu/h
614  Vmath::Vdiv(nq,physin[1],1,physin[0],1,physout[1],1);
615 
616  // v = hv/ v
617  Vmath::Vdiv(nq,physin[2],1,physin[0],1,physout[2],1);
618  }
619  }
620 
621 
622  void LinearSWE::v_ConservativeToPrimitive( )
623  {
624  int nq = GetTotPoints();
625 
626  // u = hu/h
627  Vmath::Vdiv(nq,m_fields[1]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
628 
629  // v = hv/ v
630  Vmath::Vdiv(nq,m_fields[2]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
631 
632  // \eta = h - d
633  Vmath::Vsub(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
634  }
635 
636  void LinearSWE::PrimitiveToConservative(const Array<OneD, const Array<OneD, NekDouble> >&physin,
637  Array<OneD, Array<OneD, NekDouble> >&physout)
638  {
639 
640  int nq = GetTotPoints();
641 
642  if(physin.get() == physout.get())
643  {
644  // copy indata and work with tmp array
645  Array<OneD, Array<OneD, NekDouble> >tmp(3);
646  for (int i = 0; i < 3; ++i)
647  {
648  // deep copy
649  tmp[i] = Array<OneD, NekDouble>(nq);
650  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
651  }
652 
653  // h = \eta + d
654  Vmath::Vadd(nq,tmp[0],1,m_depth,1,physout[0],1);
655 
656  // hu = h * u
657  Vmath::Vmul(nq,physout[0],1,tmp[1],1,physout[1],1);
658 
659  // hv = h * v
660  Vmath::Vmul(nq,physout[0],1,tmp[2],1,physout[2],1);
661 
662  }
663  else
664  {
665  // h = \eta + d
666  Vmath::Vadd(nq,physin[0],1,m_depth,1,physout[0],1);
667 
668  // hu = h * u
669  Vmath::Vmul(nq,physout[0],1,physin[1],1,physout[1],1);
670 
671  // hv = h * v
672  Vmath::Vmul(nq,physout[0],1,physin[2],1,physout[2],1);
673 
674  }
675 
676  }
677 
678  void LinearSWE::v_PrimitiveToConservative( )
679  {
680  int nq = GetTotPoints();
681 
682  // h = \eta + d
683  Vmath::Vadd(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
684 
685  // hu = h * u
686  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[1]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
687 
688  // hv = h * v
689  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[2]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
690  }
691 
692 
693  /**
694  * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
695  * \f$ h\mathbf{v} \f$.
696  *
697  * @param physfield Velocity field.
698  * @param velocity Velocity field.
699  */
700  void LinearSWE::GetVelocityVector(
701  const Array<OneD, Array<OneD, NekDouble> > &physfield,
702  Array<OneD, Array<OneD, NekDouble> > &velocity)
703  {
704  const int npts = physfield[0].num_elements();
705 
706  for (int i = 0; i < m_spacedim; ++i)
707  {
708  Vmath::Vcopy(npts, physfield[1+i], 1, velocity[i], 1);
709  }
710  }
711 
712 
713  void LinearSWE::v_GenerateSummary(SolverUtils::SummaryList& s)
714  {
715  ShallowWaterSystem::v_GenerateSummary(s);
716  if (m_session->DefinesSolverInfo("UpwindType"))
717  {
718  std::string UpwindType;
719  UpwindType = m_session->GetSolverInfo("UpwindType");
720  if (UpwindType == "LinearAverage")
721  {
722  SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear Average");
723  }
724  if (UpwindType == "LinearHLL")
725  {
726  SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear HLL");
727  }
728  }
729  SolverUtils::AddSummaryItem(s, "Variables", "eta should be in field[0]");
730  SolverUtils::AddSummaryItem(s, "", "u should be in field[1]");
731  SolverUtils::AddSummaryItem(s, "", "v should be in field[2]");
732  }
733 
734 } //end of namespace
735 
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2