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,
277  outarray, time);
278  //-------------------------------------------------------
279 
280 
281  //-------------------------------------------------------
282  // negate the outarray since moving terms to the rhs
283  for(i = 0; i < nvariables; ++i)
284  {
285  Vmath::Neg(nq,outarray[i],1);
286  }
287  //-------------------------------------------------------
288 
289 
290  //-------------------------------------------------
291  // Add "source terms"
292  // Input and output in physical space
293 
294  // Coriolis forcing
295  if (m_coriolis.num_elements() != 0)
296  {
297  AddCoriolis(inarray,outarray);
298  }
299 
300  // Variable Depth
301  if (m_constantDepth != true)
302  {
303  AddVariableDepth(inarray,outarray);
304  }
305  //-------------------------------------------------
306 
307  }
308  break;
309  case MultiRegions::eGalerkin:
310  case MultiRegions::eMixed_CG_Discontinuous:
311  {
312 
313  //-------------------------------------------------------
314  // Compute the fluxvector in physical space
315  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
316  fluxvector(nvariables);
317 
318  for (i = 0; i < nvariables; ++i)
319  {
320  fluxvector[i] = Array<OneD, Array<OneD, NekDouble> >(ndim);
321  for(j = 0; j < ndim; ++j)
322  {
323  fluxvector[i][j] = Array<OneD, NekDouble>(nq);
324  }
325  }
326 
327  NonlinearSWE::GetFluxVector(inarray, fluxvector);
328  //-------------------------------------------------------
329 
330 
331 
332  //-------------------------------------------------------
333  // Take the derivative of the flux terms
334  // and negate the outarray since moving terms to the rhs
335  Array<OneD,NekDouble> tmp(nq);
336  Array<OneD, NekDouble>tmp1(nq);
337 
338  for(i = 0; i < nvariables; ++i)
339  {
340  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[0],fluxvector[i][0],tmp);
341  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[1],fluxvector[i][1],tmp1);
342  Vmath::Vadd(nq,tmp,1,tmp1,1,outarray[i],1);
343  Vmath::Neg(nq,outarray[i],1);
344  }
345 
346 
347  //-------------------------------------------------
348  // Add "source terms"
349  // Input and output in physical space
350 
351  // Coriolis forcing
352  if (m_coriolis.num_elements() != 0)
353  {
354  AddCoriolis(inarray,outarray);
355  }
356 
357  // Variable Depth
358  if (m_constantDepth != true)
359  {
360  AddVariableDepth(inarray,outarray);
361  }
362  //-------------------------------------------------
363  }
364  break;
365  default:
366  ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
367  break;
368  }
369  }
370 
371 
372  void NonlinearSWE::DoOdeProjection(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
373  Array<OneD, Array<OneD, NekDouble> >&outarray,
374  const NekDouble time)
375  {
376  int i;
377  int nvariables = inarray.num_elements();
378 
379 
380  switch(m_projectionType)
381  {
382  case MultiRegions::eDiscontinuous:
383  {
384 
385  // Just copy over array
386  int npoints = GetNpoints();
387 
388  for(i = 0; i < nvariables; ++i)
389  {
390  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
391  }
392  SetBoundaryConditions(outarray, time);
393  break;
394  }
395  case MultiRegions::eGalerkin:
396  case MultiRegions::eMixed_CG_Discontinuous:
397  {
398 
399  EquationSystem::SetBoundaryConditions(time);
400  Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
401 
402  for(i = 0; i < nvariables; ++i)
403  {
404  m_fields[i]->FwdTrans(inarray[i],coeffs);
405  m_fields[i]->BwdTrans_IterPerExp(coeffs,outarray[i]);
406  }
407  break;
408  }
409  default:
410  ASSERTL0(false,"Unknown projection scheme");
411  break;
412  }
413  }
414 
415 
416  //----------------------------------------------------
417  void NonlinearSWE::SetBoundaryConditions(
418  Array<OneD, Array<OneD, NekDouble> > &inarray,
419  NekDouble time)
420  {
421  std::string varName;
422  int nvariables = m_fields.num_elements();
423  int cnt = 0;
424 
425  // Loop over Boundary Regions
426  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
427  {
428 
429  // Wall Boundary Condition
430  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
431  SpatialDomains::eWall)
432  {
433  WallBoundary2D(n, cnt, inarray);
434  }
435 
436  // Time Dependent Boundary Condition (specified in meshfile)
437  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
438  SpatialDomains::eTimeDependent)
439  {
440  for (int i = 0; i < nvariables; ++i)
441  {
442  varName = m_session->GetVariable(i);
443  m_fields[i]->EvaluateBoundaryConditions(time, varName);
444  }
445  }
446  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
447  }
448  }
449 
450  //----------------------------------------------------
451  /**
452  * @brief Wall boundary condition.
453  */
454  void NonlinearSWE::WallBoundary(
455  int bcRegion,
456  int cnt,
457  Array<OneD, Array<OneD, NekDouble> > &physarray)
458  {
459  int i;
460  int nTracePts = GetTraceTotPoints();
461  int nvariables = physarray.num_elements();
462 
463  // get physical values of the forward trace
464  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
465  for (i = 0; i < nvariables; ++i)
466  {
467  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
468  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
469  }
470 
471  // Adjust the physical values of the trace to take
472  // user defined boundaries into account
473  int e, id1, id2, npts;
474 
475  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]
476  ->GetExpSize(); ++e)
477  {
478  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
479  GetExp(e)->GetTotPoints();
480  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
481  GetPhys_Offset(e);
482  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
483  m_fields[0]->GetTraceMap()->
484  GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
485 
486  // For 2D/3D, define: v* = v - 2(v.n)n
487  Array<OneD, NekDouble> tmp(npts, 0.0);
488 
489  // Calculate (v.n)
490  for (i = 0; i < m_spacedim; ++i)
491  {
492  Vmath::Vvtvp(npts,
493  &Fwd[1+i][id2], 1,
494  &m_traceNormals[i][id2], 1,
495  &tmp[0], 1,
496  &tmp[0], 1);
497  }
498 
499  // Calculate 2.0(v.n)
500  Vmath::Smul(npts, -2.0, &tmp[0], 1, &tmp[0], 1);
501 
502  // Calculate v* = v - 2.0(v.n)n
503  for (i = 0; i < m_spacedim; ++i)
504  {
505  Vmath::Vvtvp(npts,
506  &tmp[0], 1,
507  &m_traceNormals[i][id2], 1,
508  &Fwd[1+i][id2], 1,
509  &Fwd[1+i][id2], 1);
510  }
511 
512  // copy boundary adjusted values into the boundary expansion
513  for (i = 0; i < nvariables; ++i)
514  {
515  Vmath::Vcopy(npts, &Fwd[i][id2], 1,
516  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
517  UpdatePhys())[id1], 1);
518  }
519  }
520  }
521 
522 
523  void NonlinearSWE::WallBoundary2D(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &physarray)
524  {
525 
526  int i;
527  int nTraceNumPoints = GetTraceTotPoints();
528  int nvariables = physarray.num_elements();
529 
530  // get physical values of the forward trace
531  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
532  for (i = 0; i < nvariables; ++i)
533  {
534  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
535  m_fields[i]->ExtractTracePhys(physarray[i],Fwd[i]);
536  }
537 
538  // Adjust the physical values of the trace to take
539  // user defined boundaries into account
540  int e, id1, id2, npts;
541 
542  for(e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize(); ++e)
543  {
544  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExp(e)->GetNumPoints(0);
545  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e) ;
546  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
547 
548  switch(m_expdim)
549  {
550  case 1:
551  {
552  // negate the forward flux
553  Vmath::Neg(npts,&Fwd[1][id2],1);
554  }
555  break;
556  case 2:
557  {
558  Array<OneD, NekDouble> tmp_n(npts);
559  Array<OneD, NekDouble> tmp_t(npts);
560 
561  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[0][id2],1,&tmp_n[0],1);
562  Vmath::Vvtvp(npts,&Fwd[2][id2],1,&m_traceNormals[1][id2],1,&tmp_n[0],1,&tmp_n[0],1);
563 
564  Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[1][id2],1,&tmp_t[0],1);
565  Vmath::Vvtvm(npts,&Fwd[2][id2],1,&m_traceNormals[0][id2],1,&tmp_t[0],1,&tmp_t[0],1);
566 
567  // negate the normal flux
568  Vmath::Neg(npts,tmp_n,1);
569 
570  // rotate back to Cartesian
571  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[1][id2],1,&Fwd[1][id2],1);
572  Vmath::Vvtvm(npts,&tmp_n[0],1,&m_traceNormals[0][id2],1,&Fwd[1][id2],1,&Fwd[1][id2],1);
573 
574  Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[0][id2],1,&Fwd[2][id2],1);
575  Vmath::Vvtvp(npts,&tmp_n[0],1,&m_traceNormals[1][id2],1,&Fwd[2][id2],1,&Fwd[2][id2],1);
576  }
577  break;
578  case 3:
579  ASSERTL0(false,"3D not implemented for Shallow Water Equations");
580  break;
581  default:
582  ASSERTL0(false,"Illegal expansion dimension");
583  }
584 
585 
586 
587  // copy boundary adjusted values into the boundary expansion
588  for (i = 0; i < nvariables; ++i)
589  {
590  Vmath::Vcopy(npts,&Fwd[i][id2], 1,&(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1],1);
591  }
592  }
593  }
594 
595 
596  // Physfield in conservative Form
597  void NonlinearSWE::GetFluxVector(
598  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
599  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
600  {
601  int i, j;
602  int nq = m_fields[0]->GetTotPoints();
603 
604  NekDouble g = m_g;
605  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
606 
607  // Flux vector for the mass equation
608  for (i = 0; i < m_spacedim; ++i)
609  {
610  velocity[i] = Array<OneD, NekDouble>(nq);
611  Vmath::Vcopy(nq, physfield[i+1], 1, flux[0][i], 1);
612  }
613 
614  GetVelocityVector(physfield, velocity);
615 
616  // Put (0.5 g h h) in tmp
617  Array<OneD, NekDouble> tmp(nq);
618  Vmath::Vmul(nq, physfield[0], 1, physfield[0], 1, tmp, 1);
619  Vmath::Smul(nq, 0.5*g, tmp, 1, tmp, 1);
620 
621  // Flux vector for the momentum equations
622  for (i = 0; i < m_spacedim; ++i)
623  {
624  for (j = 0; j < m_spacedim; ++j)
625  {
626  Vmath::Vmul(nq, velocity[j], 1, physfield[i+1], 1,
627  flux[i+1][j], 1);
628  }
629 
630  // Add (0.5 g h h) to appropriate field
631  Vmath::Vadd(nq, flux[i+1][i], 1, tmp, 1, flux[i+1][i], 1);
632  }
633 
634  }
635 
636  void NonlinearSWE::ConservativeToPrimitive(const Array<OneD, const Array<OneD, NekDouble> >&physin,
637  Array<OneD, Array<OneD, NekDouble> >&physout)
638  {
639  int nq = GetTotPoints();
640 
641  if(physin.get() == physout.get())
642  {
643  // copy indata and work with tmp array
644  Array<OneD, Array<OneD, NekDouble> >tmp(3);
645  for (int i = 0; i < 3; ++i)
646  {
647  // deep copy
648  tmp[i] = Array<OneD, NekDouble>(nq);
649  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
650  }
651 
652  // \eta = h - d
653  Vmath::Vsub(nq,tmp[0],1,m_depth,1,physout[0],1);
654 
655  // u = hu/h
656  Vmath::Vdiv(nq,tmp[1],1,tmp[0],1,physout[1],1);
657 
658  // v = hv/ v
659  Vmath::Vdiv(nq,tmp[2],1,tmp[0],1,physout[2],1);
660  }
661  else
662  {
663  // \eta = h - d
664  Vmath::Vsub(nq,physin[0],1,m_depth,1,physout[0],1);
665 
666  // u = hu/h
667  Vmath::Vdiv(nq,physin[1],1,physin[0],1,physout[1],1);
668 
669  // v = hv/ v
670  Vmath::Vdiv(nq,physin[2],1,physin[0],1,physout[2],1);
671  }
672  }
673 
674 
675  void NonlinearSWE::v_ConservativeToPrimitive( )
676  {
677  int nq = GetTotPoints();
678 
679  // u = hu/h
680  Vmath::Vdiv(nq,m_fields[1]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
681 
682  // v = hv/ v
683  Vmath::Vdiv(nq,m_fields[2]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
684 
685  // \eta = h - d
686  Vmath::Vsub(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
687  }
688 
689  void NonlinearSWE::PrimitiveToConservative(const Array<OneD, const Array<OneD, NekDouble> >&physin,
690  Array<OneD, Array<OneD, NekDouble> >&physout)
691  {
692 
693  int nq = GetTotPoints();
694 
695  if(physin.get() == physout.get())
696  {
697  // copy indata and work with tmp array
698  Array<OneD, Array<OneD, NekDouble> >tmp(3);
699  for (int i = 0; i < 3; ++i)
700  {
701  // deep copy
702  tmp[i] = Array<OneD, NekDouble>(nq);
703  Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
704  }
705 
706  // h = \eta + d
707  Vmath::Vadd(nq,tmp[0],1,m_depth,1,physout[0],1);
708 
709  // hu = h * u
710  Vmath::Vmul(nq,physout[0],1,tmp[1],1,physout[1],1);
711 
712  // hv = h * v
713  Vmath::Vmul(nq,physout[0],1,tmp[2],1,physout[2],1);
714 
715  }
716  else
717  {
718  // h = \eta + d
719  Vmath::Vadd(nq,physin[0],1,m_depth,1,physout[0],1);
720 
721  // hu = h * u
722  Vmath::Vmul(nq,physout[0],1,physin[1],1,physout[1],1);
723 
724  // hv = h * v
725  Vmath::Vmul(nq,physout[0],1,physin[2],1,physout[2],1);
726 
727  }
728 
729  }
730 
731  void NonlinearSWE::v_PrimitiveToConservative( )
732  {
733  int nq = GetTotPoints();
734 
735  // h = \eta + d
736  Vmath::Vadd(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
737 
738  // hu = h * u
739  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[1]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
740 
741  // hv = h * v
742  Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[2]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
743  }
744 
745 
746  /**
747  * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
748  * \f$ h\mathbf{v} \f$.
749  *
750  * @param physfield Momentum field.
751  * @param velocity Velocity field.
752  */
753  void NonlinearSWE::GetVelocityVector(
754  const Array<OneD, Array<OneD, NekDouble> > &physfield,
755  Array<OneD, Array<OneD, NekDouble> > &velocity)
756  {
757  const int npts = physfield[0].num_elements();
758 
759  for (int i = 0; i < m_spacedim; ++i)
760  {
761  Vmath::Vdiv(npts, physfield[1+i], 1, physfield[0], 1,
762  velocity[i], 1);
763  }
764  }
765 
766 
767  void NonlinearSWE::v_GenerateSummary(SolverUtils::SummaryList& s)
768  {
769  ShallowWaterSystem::v_GenerateSummary(s);
770  SolverUtils::AddSummaryItem(s, "Variables", "h should be in field[0]");
771  SolverUtils::AddSummaryItem(s, "", "hu should be in field[1]");
772  SolverUtils::AddSummaryItem(s, "", "hv should be in field[2]");
773  }
774 
775 } //end of namespace
776 
Generated on Wed Mar 4 2015 16:33:16 for Nektar++ by   doxygen 1.8.1.2