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