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,
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;
271  case MultiRegions::eGalerkin:
272  case MultiRegions::eMixed_CG_Discontinuous:
273  {
274 
275  //-------------------------------------------------------
276  // Compute the fluxvector in physical space
277  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
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 
326  void LinearSWE::DoOdeProjection(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
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  {
336  case MultiRegions::eDiscontinuous:
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  }
349  case MultiRegions::eGalerkin:
350  case MultiRegions::eMixed_CG_Discontinuous:
351  {
352 
353  EquationSystem::SetBoundaryConditions(time);
354  Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
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  //----------------------------------------------------
371  void LinearSWE::SetBoundaryConditions(
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 (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
384  SpatialDomains::eWall)
385  {
386  WallBoundary2D(n, cnt, inarray);
387  }
388 
389  // Time Dependent Boundary Condition (specified in meshfile)
390  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
391  SpatialDomains::eTimeDependent)
392  {
393  for (int i = 0; i < nvariables; ++i)
394  {
395  varName = m_session->GetVariable(i);
396  m_fields[i]->EvaluateBoundaryConditions(time, varName);
397  }
398  }
399  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
400  }
401  }
402 
403  //----------------------------------------------------
404  /**
405  * @brief Wall boundary condition.
406  */
407  void LinearSWE::WallBoundary(
408  int bcRegion,
409  int cnt,
410  Array<OneD, Array<OneD, NekDouble> > &physarray)
411  {
412  int i;
413  int nTracePts = GetTraceTotPoints();
414  int nvariables = physarray.num_elements();
415 
416  // get physical values of the forward trace
417  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
418  for (i = 0; i < nvariables; ++i)
419  {
420  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
421  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
422  }
423 
424  // Adjust the physical values of the trace to take
425  // user defined boundaries into account
426  int e, id1, id2, npts;
427 
428  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]
429  ->GetExpSize(); ++e)
430  {
431  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
432  GetExp(e)->GetTotPoints();
433  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
434  GetPhys_Offset(e);
435  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
436  m_fields[0]->GetTraceMap()->
437  GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
438 
439  // For 2D/3D, define: v* = v - 2(v.n)n
440  Array<OneD, NekDouble> tmp(npts, 0.0);
441 
442  // Calculate (v.n)
443  for (i = 0; i < m_spacedim; ++i)
444  {
445  Vmath::Vvtvp(npts,
446  &Fwd[1+i][id2], 1,
447  &m_traceNormals[i][id2], 1,
448  &tmp[0], 1,
449  &tmp[0], 1);
450  }
451 
452  // Calculate 2.0(v.n)
453  Vmath::Smul(npts, -2.0, &tmp[0], 1, &tmp[0], 1);
454 
455  // Calculate v* = v - 2.0(v.n)n
456  for (i = 0; i < m_spacedim; ++i)
457  {
458  Vmath::Vvtvp(npts,
459  &tmp[0], 1,
460  &m_traceNormals[i][id2], 1,
461  &Fwd[1+i][id2], 1,
462  &Fwd[1+i][id2], 1);
463  }
464 
465  // copy boundary adjusted values into the boundary expansion
466  for (i = 0; i < nvariables; ++i)
467  {
468  Vmath::Vcopy(npts, &Fwd[i][id2], 1,
469  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
470  UpdatePhys())[id1], 1);
471  }
472  }
473  }
474 
475 
476  void LinearSWE::WallBoundary2D(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &physarray)
477  {
478 
479  int i;
480  int nTraceNumPoints = GetTraceTotPoints();
481  int nvariables = physarray.num_elements();
482 
483  // get physical values of the forward trace
484  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
485  for (i = 0; i < nvariables; ++i)
486  {
487  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
488  m_fields[i]->ExtractTracePhys(physarray[i],Fwd[i]);
489  }
490 
491  // Adjust the physical values of the trace to take
492  // user defined boundaries into account
493  int e, id1, id2, npts;
494 
495  for(e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize(); ++e)
496  {
497  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExp(e)->GetNumPoints(0);
498  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e) ;
499  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
500 
501  switch(m_expdim)
502  {
503  case 1:
504  {
505  // negate the forward flux
506  Vmath::Neg(npts,&Fwd[1][id2],1);
507  }
508  break;
509  case 2:
510  {
511  Array<OneD, NekDouble> tmp_n(npts);
512  Array<OneD, NekDouble> tmp_t(npts);
513 
514  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[0][id2],1,&tmp_n[0],1);
515  Vmath::Vvtvp(npts,&Fwd[2][id2],1,&m_traceNormals[1][id2],1,&tmp_n[0],1,&tmp_n[0],1);
516 
517  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[1][id2],1,&tmp_t[0],1);
518  Vmath::Vvtvm(npts,&Fwd[2][id2],1,&m_traceNormals[0][id2],1,&tmp_t[0],1,&tmp_t[0],1);
519 
520  // negate the normal flux
521  Vmath::Neg(npts,tmp_n,1);
522 
523  // rotate back to Cartesian
524  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[1][id2],1,&Fwd[1][id2],1);
525  Vmath::Vvtvm(npts,&tmp_n[0],1,&m_traceNormals[0][id2],1,&Fwd[1][id2],1,&Fwd[1][id2],1);
526 
527  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[0][id2],1,&Fwd[2][id2],1);
528  Vmath::Vvtvp(npts,&tmp_n[0],1,&m_traceNormals[1][id2],1,&Fwd[2][id2],1,&Fwd[2][id2],1);
529  }
530  break;
531  case 3:
532  ASSERTL0(false,"3D not implemented for Shallow Water Equations");
533  break;
534  default:
535  ASSERTL0(false,"Illegal expansion dimension");
536  }
537 
538 
539 
540  // copy boundary adjusted values into the boundary expansion
541  for (i = 0; i < nvariables; ++i)
542  {
543  Vmath::Vcopy(npts,&Fwd[i][id2], 1,&(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1],1);
544  }
545  }
546  }
547 
548 
549  // Physfield in primitive Form
550  void LinearSWE::GetFluxVector(
551  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
552  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
553  {
554  int i, j;
555  int nq = m_fields[0]->GetTotPoints();
556 
557  NekDouble g = m_g;
558 
559  // Flux vector for the mass equation
560  for (i = 0; i < m_spacedim; ++i)
561  {
562  Vmath::Vmul(nq, m_depth, 1, physfield[i+1], 1, flux[0][i], 1);
563  }
564 
565  // Put (g eta) in tmp
566  Array<OneD, NekDouble> tmp(nq);
567  Vmath::Smul(nq, g, physfield[0], 1, tmp, 1);
568 
569  // Flux vector for the momentum equations
570  for (i = 0; i < m_spacedim; ++i)
571  {
572  for (j = 0; j < m_spacedim; ++j)
573  {
574  // must zero fluxes as not initialised to zero in AdvectionWeakDG ...
575  Vmath::Zero(nq, flux[i+1][j], 1);
576  }
577 
578  // Add (g eta) to appropriate field
579  Vmath::Vadd(nq, flux[i+1][i], 1, tmp, 1, flux[i+1][i], 1);
580  }
581 
582  }
583 
584  void LinearSWE::ConservativeToPrimitive(const Array<OneD, const Array<OneD, NekDouble> >&physin,
585  Array<OneD, Array<OneD, NekDouble> >&physout)
586  {
587  int nq = GetTotPoints();
588 
589  if(physin.get() == physout.get())
590  {
591  // copy indata and work with tmp array
592  Array<OneD, Array<OneD, NekDouble> >tmp(3);
593  for (int i = 0; i < 3; ++i)
594  {
595  // deep copy
596  tmp[i] = Array<OneD, NekDouble>(nq);
597  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
598  }
599 
600  // \eta = h - d
601  Vmath::Vsub(nq,tmp[0],1,m_depth,1,physout[0],1);
602 
603  // u = hu/h
604  Vmath::Vdiv(nq,tmp[1],1,tmp[0],1,physout[1],1);
605 
606  // v = hv/ v
607  Vmath::Vdiv(nq,tmp[2],1,tmp[0],1,physout[2],1);
608  }
609  else
610  {
611  // \eta = h - d
612  Vmath::Vsub(nq,physin[0],1,m_depth,1,physout[0],1);
613 
614  // u = hu/h
615  Vmath::Vdiv(nq,physin[1],1,physin[0],1,physout[1],1);
616 
617  // v = hv/ v
618  Vmath::Vdiv(nq,physin[2],1,physin[0],1,physout[2],1);
619  }
620  }
621 
622 
623  void LinearSWE::v_ConservativeToPrimitive( )
624  {
625  int nq = GetTotPoints();
626 
627  // u = hu/h
628  Vmath::Vdiv(nq,m_fields[1]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
629 
630  // v = hv/ v
631  Vmath::Vdiv(nq,m_fields[2]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
632 
633  // \eta = h - d
634  Vmath::Vsub(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
635  }
636 
637  void LinearSWE::PrimitiveToConservative(const Array<OneD, const Array<OneD, NekDouble> >&physin,
638  Array<OneD, Array<OneD, NekDouble> >&physout)
639  {
640 
641  int nq = GetTotPoints();
642 
643  if(physin.get() == physout.get())
644  {
645  // copy indata and work with tmp array
646  Array<OneD, Array<OneD, NekDouble> >tmp(3);
647  for (int i = 0; i < 3; ++i)
648  {
649  // deep copy
650  tmp[i] = Array<OneD, NekDouble>(nq);
651  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
652  }
653 
654  // h = \eta + d
655  Vmath::Vadd(nq,tmp[0],1,m_depth,1,physout[0],1);
656 
657  // hu = h * u
658  Vmath::Vmul(nq,physout[0],1,tmp[1],1,physout[1],1);
659 
660  // hv = h * v
661  Vmath::Vmul(nq,physout[0],1,tmp[2],1,physout[2],1);
662 
663  }
664  else
665  {
666  // h = \eta + d
667  Vmath::Vadd(nq,physin[0],1,m_depth,1,physout[0],1);
668 
669  // hu = h * u
670  Vmath::Vmul(nq,physout[0],1,physin[1],1,physout[1],1);
671 
672  // hv = h * v
673  Vmath::Vmul(nq,physout[0],1,physin[2],1,physout[2],1);
674 
675  }
676 
677  }
678 
679  void LinearSWE::v_PrimitiveToConservative( )
680  {
681  int nq = GetTotPoints();
682 
683  // h = \eta + d
684  Vmath::Vadd(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
685 
686  // hu = h * u
687  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[1]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
688 
689  // hv = h * v
690  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[2]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
691  }
692 
693 
694  /**
695  * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
696  * \f$ h\mathbf{v} \f$.
697  *
698  * @param physfield Velocity field.
699  * @param velocity Velocity field.
700  */
701  void LinearSWE::GetVelocityVector(
702  const Array<OneD, Array<OneD, NekDouble> > &physfield,
703  Array<OneD, Array<OneD, NekDouble> > &velocity)
704  {
705  const int npts = physfield[0].num_elements();
706 
707  for (int i = 0; i < m_spacedim; ++i)
708  {
709  Vmath::Vcopy(npts, physfield[1+i], 1, velocity[i], 1);
710  }
711  }
712 
713 
714  void LinearSWE::v_GenerateSummary(SolverUtils::SummaryList& s)
715  {
716  ShallowWaterSystem::v_GenerateSummary(s);
717  if (m_session->DefinesSolverInfo("UpwindType"))
718  {
719  std::string UpwindType;
720  UpwindType = m_session->GetSolverInfo("UpwindType");
721  if (UpwindType == "LinearAverage")
722  {
723  SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear Average");
724  }
725  if (UpwindType == "LinearHLL")
726  {
727  SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear HLL");
728  }
729  }
730  SolverUtils::AddSummaryItem(s, "Variables", "eta should be in field[0]");
731  SolverUtils::AddSummaryItem(s, "", "u should be in field[1]");
732  SolverUtils::AddSummaryItem(s, "", "v should be in field[2]");
733  }
734 
735 } //end of namespace
736 
Generated on Wed Mar 4 2015 16:33:15 for Nektar++ by   doxygen 1.8.1.2