LinearSWE.cpp 22.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
///////////////////////////////////////////////////////////////////////////////
//
// File LinearSWE.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
Claes Eskilsson's avatar
Claes Eskilsson committed
32
// Description: Linear Shallow water equations in primitive variables
33 34 35
//
///////////////////////////////////////////////////////////////////////////////

Claes Eskilsson's avatar
Claes Eskilsson committed
36 37 38 39
#include <iostream>
#include <iomanip>
#include <boost/algorithm/string.hpp>

40
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
41 42 43 44
#include <ShallowWaterSolver/EquationSystems/LinearSWE.h>

namespace Nektar
{
Claes Eskilsson's avatar
Claes Eskilsson committed
45 46 47 48
  string LinearSWE::className = 
     SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
	  "LinearSWE", LinearSWE::create, 
          "Linear shallow water equation in primitive variables.");
49
  
Claes Eskilsson's avatar
Claes Eskilsson committed
50 51
  LinearSWE::LinearSWE(
          const LibUtilities::SessionReaderSharedPtr& pSession)
52
    : ShallowWaterSystem(pSession)
53
  {
54 55 56 57 58
  }

  void LinearSWE::v_InitObject()
  {
      ShallowWaterSystem::v_InitObject();
59 60 61 62 63 64 65 66 67 68 69
      
    if (m_explicitAdvection)
      {
	m_ode.DefineOdeRhs     (&LinearSWE::DoOdeRhs,        this);
	m_ode.DefineProjection (&LinearSWE::DoOdeProjection, this);
      }
    else
      {
	ASSERTL0(false, "Implicit SWE not set up.");
      }

Claes Eskilsson's avatar
Claes Eskilsson committed
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
    // Type of advection class to be used
    switch(m_projectionType)
      {
	// Continuous field 
      case MultiRegions::eGalerkin:
	{
	  //  Do nothing 
	  break;
	}
	// Discontinuous field 
      case MultiRegions::eDiscontinuous:
	{
	  string advName;
	  string diffName;
	  string riemName;
          
	  //---------------------------------------------------------------
	  // Setting up advection and diffusion operators
	  // NB: diffusion not set up for SWE at the moment
	  //     but kept here for future use ...
	  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
	  // m_session->LoadSolverInfo("DiffusionType", diffName, "LDGEddy");
	  m_advection = SolverUtils::GetAdvectionFactory()
	    .CreateInstance(advName, advName);
	  // m_diffusion = SolverUtils::GetDiffusionFactory()
	  //                             .CreateInstance(diffName, diffName);
	  
	  m_advection->SetFluxVector(&LinearSWE::GetFluxVector, this);
	  // m_diffusion->SetFluxVectorNS(&ShallowWaterSystem::
	  //                                  GetEddyViscosityFluxVector, this); 
         
	  // Setting up Riemann solver for advection operator
	  m_session->LoadSolverInfo("UpwindType", riemName, "NoSolver");
	  if ((riemName == "LinearHLL") && (m_constantDepth != true))
	    {
	      ASSERTL0(false,"LinearHLL only valid for constant depth"); 
	    }
	  m_riemannSolver = SolverUtils::GetRiemannSolverFactory()
	    .CreateInstance(riemName);
         
       	  // Setting up upwind solver for diffusion operator
	  // m_riemannSolverLDG = SolverUtils::GetRiemannSolverFactory()
	  //                                 .CreateInstance("UpwindLDG");
	  
	  // Setting up parameters for advection operator Riemann solver 
Dave Moxey's avatar
Dave Moxey committed
115
	  m_riemannSolver->SetParam (
Claes Eskilsson's avatar
Claes Eskilsson committed
116 117
                                     "gravity",  
                                     &LinearSWE::GetGravity,   this);
118
	  m_riemannSolver->SetAuxiliary(
Claes Eskilsson's avatar
Claes Eskilsson committed
119 120
                                     "velLoc", 
                                     &LinearSWE::GetVelLoc,  this);
Dave Moxey's avatar
Dave Moxey committed
121
	  m_riemannSolver->SetVector(
Claes Eskilsson's avatar
Claes Eskilsson committed
122 123 124 125 126 127 128 129 130
				     "N",
				     &LinearSWE::GetNormals, this);
	
	  // The numerical flux for linear SWE requires depth information
	  int nTracePointsTot = m_fields[0]->GetTrace()->GetTotPoints();
	  m_dFwd = Array<OneD, NekDouble>(nTracePointsTot);
	  m_dBwd = Array<OneD, NekDouble>(nTracePointsTot);
	  m_fields[0]->GetFwdBwdTracePhys(m_depth, m_dFwd, m_dBwd);
	  CopyBoundaryTrace(m_dFwd,m_dBwd);
Dave Moxey's avatar
Dave Moxey committed
131
	  m_riemannSolver->SetScalar(
Claes Eskilsson's avatar
Claes Eskilsson committed
132 133
				     "depthFwd",
				     &LinearSWE::GetDepthFwd, this);
Dave Moxey's avatar
Dave Moxey committed
134
	  m_riemannSolver->SetScalar(
Claes Eskilsson's avatar
Claes Eskilsson committed
135 136 137 138 139 140 141
				     "depthBwd",
				     &LinearSWE::GetDepthBwd, this);

	  // Setting up parameters for diffusion operator Riemann solver
	  // m_riemannSolverLDG->AddParam (
	  //                     "gravity",  
	  //                     &NonlinearSWE::GetGravity,   this);
142
	  // m_riemannSolverLDG->AddAuxiliary(
Claes Eskilsson's avatar
Claes Eskilsson committed
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
	  //                     "velLoc", 
	  //                     &NonlinearSWE::GetVelLoc,  this);
	  // m_riemannSolverLDG->AddVector(
	  //                     "N",
	  //                     &NonlinearSWE::GetNormals, this);
          
	  // Concluding initialisation of advection / diffusion operators
	  m_advection->SetRiemannSolver   (m_riemannSolver);
	  //m_diffusion->SetRiemannSolver   (m_riemannSolverLDG);
	  m_advection->InitObject         (m_session, m_fields);
	  //m_diffusion->InitObject         (m_session, m_fields);
	  break;
	}
      default:
	{
	  ASSERTL0(false, "Unsupported projection type.");
	  break;
	}
      }
    

164 165 166 167 168 169 170
  }
  
  LinearSWE::~LinearSWE()
  {
    
  }
 
Claes Eskilsson's avatar
Claes Eskilsson committed
171
  // physarray contains the conservative variables
172
  void LinearSWE::AddCoriolis(const Array<OneD, const Array<OneD, NekDouble> > &physarray,
Claes Eskilsson's avatar
Claes Eskilsson committed
173
			      Array<OneD, Array<OneD, NekDouble> > &outarray)
174 175
  {
    
Claes Eskilsson's avatar
Claes Eskilsson committed
176 177
    int ncoeffs    = GetNcoeffs();
    int nq         = GetTotPoints();
178
    
Claes Eskilsson's avatar
Claes Eskilsson committed
179 180 181
    Array<OneD, NekDouble> tmp(nq);
    Array<OneD, NekDouble> mod(ncoeffs);
	
182 183
    switch(m_projectionType)
      {
184
      case MultiRegions::eDiscontinuous:
Claes Eskilsson's avatar
Claes Eskilsson committed
185
	{	  
186 187
	  // add to u equation
	  Vmath::Vmul(nq,m_coriolis,1,physarray[2],1,tmp,1);
Claes Eskilsson's avatar
Claes Eskilsson committed
188 189 190 191 192
	  m_fields[0]->IProductWRTBase(tmp,mod);
	  m_fields[0]->MultiplyByElmtInvMass(mod,mod);
	  m_fields[0]->BwdTrans(mod,tmp);
	  Vmath::Vadd(nq,tmp,1,outarray[1],1,outarray[1],1);
	   
193 194 195
	  // add to v equation
	  Vmath::Vmul(nq,m_coriolis,1,physarray[1],1,tmp,1);
	  Vmath::Neg(nq,tmp,1);
Claes Eskilsson's avatar
Claes Eskilsson committed
196 197 198 199
	  m_fields[0]->IProductWRTBase(tmp,mod);
	  m_fields[0]->MultiplyByElmtInvMass(mod,mod);
	  m_fields[0]->BwdTrans(mod,tmp);
	  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 206 207 208 209 210 211 212 213 214 215
	{
	  // add to u equation
	  Vmath::Vmul(nq,m_coriolis,1,physarray[2],1,tmp,1);
	  Vmath::Vadd(nq,tmp,1,outarray[1],1,outarray[1],1);
	  
	  // add to v equation
	  Vmath::Vmul(nq,m_coriolis,1,physarray[1],1,tmp,1);
	  Vmath::Neg(nq,tmp,1);
	  Vmath::Vadd(nq,tmp,1,outarray[2],1,outarray[2],1);
	}
	break;
      default:
Claes Eskilsson's avatar
Claes Eskilsson committed
216
	ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
217 218
	break;
      }
Claes Eskilsson's avatar
Claes Eskilsson committed
219 220


221
  }
Claes Eskilsson's avatar
Claes Eskilsson committed
222
  
223
  void LinearSWE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,  
Claes Eskilsson's avatar
Claes Eskilsson committed
224 225
			            Array<OneD,       Array<OneD, NekDouble> >&outarray, 
			      const NekDouble time) 
226
  {
Claes Eskilsson's avatar
Claes Eskilsson committed
227 228
    int i, j;
    int ndim       = m_spacedim;
229 230
    int nvariables = inarray.num_elements();
    int nq         = GetTotPoints();
Claes Eskilsson's avatar
Claes Eskilsson committed
231
  
232 233 234
    
    switch(m_projectionType)
      {
235
      case MultiRegions::eDiscontinuous:
236
	{
Claes Eskilsson's avatar
Claes Eskilsson committed
237
	 
238
	  //-------------------------------------------------------
Claes Eskilsson's avatar
Claes Eskilsson committed
239 240 241 242
	  // Compute the DG advection including the numerical flux
	  // by using SolverUtils/Advection 
	  // Input and output in physical space
	  Array<OneD, Array<OneD, NekDouble> > advVel;
243
	  
Claes Eskilsson's avatar
Claes Eskilsson committed
244
	  m_advection->Advect(nvariables, m_fields, advVel, inarray, outarray);
245 246 247 248 249 250 251
	  //-------------------------------------------------------
	  
	  
	  //-------------------------------------------------------
	  // negate the outarray since moving terms to the rhs
	  for(i = 0; i < nvariables; ++i)
	    {
Claes Eskilsson's avatar
Claes Eskilsson committed
252
	      Vmath::Neg(nq,outarray[i],1);
253 254 255 256 257
	    }
	  //-------------------------------------------------------
	  
	  
	  //-------------------------------------------------
Claes Eskilsson's avatar
Claes Eskilsson committed
258 259 260 261 262
	  // Add "source terms"
	  // Input and output in physical space
	  	  
	  // Coriolis forcing
	  if (m_coriolis.num_elements() != 0)
263
	    {
Claes Eskilsson's avatar
Claes Eskilsson committed
264
	      AddCoriolis(inarray,outarray);
265
	    }
Claes Eskilsson's avatar
Claes Eskilsson committed
266 267
	  //------------------------------------------------- 
	   
268 269
	}
	break;
270
      case MultiRegions::eGalerkin:
271
      case MultiRegions::eMixed_CG_Discontinuous:
272
	{
Claes Eskilsson's avatar
Claes Eskilsson committed
273 274 275 276 277 278 279 280 281 282 283 284 285 286
	
	  //-------------------------------------------------------
	  // Compute the fluxvector in physical space
	  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > 
	    fluxvector(nvariables);
	  
	  for (i = 0; i < nvariables; ++i)
            {
	      fluxvector[i] = Array<OneD, Array<OneD, NekDouble> >(ndim);
	      for(j = 0; j < ndim; ++j)
                {
		  fluxvector[i][j] = Array<OneD, NekDouble>(nq);
                }
            }
287
	  
Claes Eskilsson's avatar
Claes Eskilsson committed
288 289 290 291 292 293 294
	  LinearSWE::GetFluxVector(inarray, fluxvector);
	  //-------------------------------------------------------

	  
	  //-------------------------------------------------------
	  // Take the derivative of the flux terms
	  // and negate the outarray since moving terms to the rhs
295
	  Array<OneD,NekDouble> tmp(nq);
Claes Eskilsson's avatar
Claes Eskilsson committed
296
	  Array<OneD, NekDouble>tmp1(nq);           
297 298 299
	  
	  for(i = 0; i < nvariables; ++i)
	    {
Claes Eskilsson's avatar
Claes Eskilsson committed
300 301
	      m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[0],fluxvector[i][0],tmp);
	      m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[1],fluxvector[i][1],tmp1);
302 303 304 305 306
	      Vmath::Vadd(nq,tmp,1,tmp1,1,outarray[i],1);
	      Vmath::Neg(nq,outarray[i],1);
	    }
	  
	  //-------------------------------------------------
Claes Eskilsson's avatar
Claes Eskilsson committed
307 308 309 310 311 312 313 314 315
	  // Add "source terms"
	  // Input and output in physical space
	  
	  // Coriolis forcing
	  if (m_coriolis.num_elements() != 0)
	    {
	      AddCoriolis(inarray,outarray);
	    }
	  //------------------------------------------------- 
316 317 318
	}
	break;
      default:
Claes Eskilsson's avatar
Claes Eskilsson committed
319
	ASSERTL0(false,"Unknown projection scheme for the NonlinearSWE");
320 321 322 323 324 325
	break;
      }
  }
 

  void LinearSWE::DoOdeProjection(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Claes Eskilsson's avatar
Claes Eskilsson committed
326 327
				           Array<OneD,       Array<OneD, NekDouble> >&outarray,
				     const NekDouble time)
328 329 330 331 332 333 334
  {
    int i;
    int nvariables = inarray.num_elements();
   
    
    switch(m_projectionType)
      {
335
      case MultiRegions::eDiscontinuous:
336
	{
Claes Eskilsson's avatar
Claes Eskilsson committed
337
	  
338 339 340 341 342
	  // Just copy over array
	  int npoints = GetNpoints();
	  
	  for(i = 0; i < nvariables; ++i)
	    {
Claes Eskilsson's avatar
Claes Eskilsson committed
343
	      Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
344
	    }
Claes Eskilsson's avatar
Claes Eskilsson committed
345 346
	  SetBoundaryConditions(outarray, time);
	  break;
347
	}
348
      case MultiRegions::eGalerkin:
349
      case MultiRegions::eMixed_CG_Discontinuous:
350
	{
Claes Eskilsson's avatar
Claes Eskilsson committed
351
	  
352 353 354 355
	  EquationSystem::SetBoundaryConditions(time);
	  Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
	  
	  for(i = 0; i < nvariables; ++i)
356 357
          {
              m_fields[i]->FwdTrans(inarray[i],coeffs);
Claes Eskilsson's avatar
Claes Eskilsson committed
358
	      m_fields[i]->BwdTrans_IterPerExp(coeffs,outarray[i]);
359
          }
360 361 362 363 364 365 366 367 368 369
	  break;
	}
      default:
	ASSERTL0(false,"Unknown projection scheme");
	break;
      }
  }
  

   //----------------------------------------------------
370 371 372 373 374 375 376
  void LinearSWE::SetBoundaryConditions(
    Array<OneD, Array<OneD, NekDouble> > &inarray,
    NekDouble time)
  { 
      std::string varName;
      int nvariables = m_fields.num_elements();
      int cnt = 0;
377

378 379
      // loop over Boundary Regions
      for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
380
      {	
381 382 383 384 385 386
          // Wall Boundary Condition
          if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() == 
              SpatialDomains::eWall)
          {
              WallBoundary2D(n, cnt, inarray);
          }
387
	
388 389 390 391 392 393 394
          // Time Dependent Boundary Condition (specified in meshfile)
          if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() == 
              SpatialDomains::eTimeDependent)
          {
              for (int i = 0; i < nvariables; ++i)
              {
                  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
      }
  }
  
  //----------------------------------------------------
Claes Eskilsson's avatar
Claes Eskilsson committed
403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474
  /**
     * @brief Wall boundary condition.
     */
    void LinearSWE::WallBoundary(
        int                                   bcRegion,
        int                                   cnt, 
        Array<OneD, Array<OneD, NekDouble> > &physarray)
    { 
        int i;
        int nTracePts = GetTraceTotPoints();
        int nvariables      = physarray.num_elements();
        
        // get physical values of the forward trace
        Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
        for (i = 0; i < nvariables; ++i)
        {
            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
        }
        
        // Adjust the physical values of the trace to take 
        // user defined boundaries into account
        int e, id1, id2, npts;
        
        for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]
                 ->GetExpSize(); ++e)
        {
            npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
                GetExp(e)->GetTotPoints();
            id1  = m_fields[0]->GetBndCondExpansions()[bcRegion]->
                GetPhys_Offset(e);
            id2  = m_fields[0]->GetTrace()->GetPhys_Offset(
                        m_fields[0]->GetTraceMap()->
                                    GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
            
            // For 2D/3D, define: v* = v - 2(v.n)n
            Array<OneD, NekDouble> tmp(npts, 0.0);

            // Calculate (v.n)
            for (i = 0; i < m_spacedim; ++i)
            {
                Vmath::Vvtvp(npts,
                             &Fwd[1+i][id2], 1,
                             &m_traceNormals[i][id2], 1,
                             &tmp[0], 1,
                             &tmp[0], 1);    
            }

            // Calculate 2.0(v.n)
            Vmath::Smul(npts, -2.0, &tmp[0], 1, &tmp[0], 1);
            
            // Calculate v* = v - 2.0(v.n)n
            for (i = 0; i < m_spacedim; ++i)
            {
                Vmath::Vvtvp(npts,
                             &tmp[0], 1,
                             &m_traceNormals[i][id2], 1,
                             &Fwd[1+i][id2], 1,
                             &Fwd[1+i][id2], 1);
            }
            
            // copy boundary adjusted values into the boundary expansion
            for (i = 0; i < nvariables; ++i)
            {
                Vmath::Vcopy(npts, &Fwd[i][id2], 1,
                             &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
                             UpdatePhys())[id1], 1);
            }
        }
    }
    

475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536
  void LinearSWE::WallBoundary2D(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &physarray)
  { 

    int i;
    int nTraceNumPoints = GetTraceTotPoints();
    int nvariables      = physarray.num_elements();
    
    // get physical values of the forward trace
    Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
    for (i = 0; i < nvariables; ++i)
      {
	Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
	m_fields[i]->ExtractTracePhys(physarray[i],Fwd[i]);
      }
    
    // Adjust the physical values of the trace to take 
    // user defined boundaries into account
    int e, id1, id2, npts;
    
    for(e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize(); ++e)
      {
	npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExp(e)->GetNumPoints(0);
	id1  = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e) ;
	id2  = m_fields[0]->GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
	
	switch(m_expdim)
	  {
	  case 1:
	    {
	      // negate the forward flux
	      Vmath::Neg(npts,&Fwd[1][id2],1);	
	    }
	    break;
	  case 2:
	    {
	      Array<OneD, NekDouble> tmp_n(npts);
	      Array<OneD, NekDouble> tmp_t(npts);
	      
	      Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[0][id2],1,&tmp_n[0],1);
	      Vmath::Vvtvp(npts,&Fwd[2][id2],1,&m_traceNormals[1][id2],1,&tmp_n[0],1,&tmp_n[0],1);
	      
	      Vmath::Vmul(npts,&Fwd[1][id2],1,&m_traceNormals[1][id2],1,&tmp_t[0],1);
	      Vmath::Vvtvm(npts,&Fwd[2][id2],1,&m_traceNormals[0][id2],1,&tmp_t[0],1,&tmp_t[0],1);
	      
	      // negate the normal flux
	      Vmath::Neg(npts,tmp_n,1);		      
	      
	      // rotate back to Cartesian
	      Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[1][id2],1,&Fwd[1][id2],1);
	      Vmath::Vvtvm(npts,&tmp_n[0],1,&m_traceNormals[0][id2],1,&Fwd[1][id2],1,&Fwd[1][id2],1);
	      
	      Vmath::Vmul(npts,&tmp_t[0],1,&m_traceNormals[0][id2],1,&Fwd[2][id2],1);
	      Vmath::Vvtvp(npts,&tmp_n[0],1,&m_traceNormals[1][id2],1,&Fwd[2][id2],1,&Fwd[2][id2],1);
	    }
	    break;
	  case 3:
	    ASSERTL0(false,"3D not implemented for Shallow Water Equations");
	    break;
	  default:
	    ASSERTL0(false,"Illegal expansion dimension");
	  }

Claes Eskilsson's avatar
Claes Eskilsson committed
537 538
	

539 540 541 542 543 544 545 546 547
	// copy boundary adjusted values into the boundary expansion
	for (i = 0; i < nvariables; ++i)
	  {
	    Vmath::Vcopy(npts,&Fwd[i][id2], 1,&(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1],1);
	  }
      }
  }
  
    
Claes Eskilsson's avatar
Claes Eskilsson committed
548 549 550 551
  // Physfield in primitive Form
 void LinearSWE::GetFluxVector(
     const Array<OneD, const Array<OneD, NekDouble> > &physfield, 
           Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
552
  {
Claes Eskilsson's avatar
Claes Eskilsson committed
553
    int i, j;
554 555 556
    int nq = m_fields[0]->GetTotPoints();
    
    NekDouble g = m_g;
Claes Eskilsson's avatar
Claes Eskilsson committed
557 558 559
   
    // Flux vector for the mass equation
    for (i = 0; i < m_spacedim; ++i)
560
      {
Claes Eskilsson's avatar
Claes Eskilsson committed
561
       	Vmath::Vmul(nq, m_depth, 1, physfield[i+1], 1, flux[0][i], 1);
562
      }
Claes Eskilsson's avatar
Claes Eskilsson committed
563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586
    
    // Put (g eta) in tmp
     Array<OneD, NekDouble> tmp(nq);
     Vmath::Smul(nq, g, physfield[0], 1, tmp, 1);
     
     // Flux vector for the momentum equations
     for (i = 0; i < m_spacedim; ++i)
       {
	 for (j = 0; j < m_spacedim; ++j)
	   {
	     // must zero fluxes as not initialised to zero in AdvectionWeakDG ... 
	     Vmath::Zero(nq, flux[i+1][j], 1);
	   }
            
	 // Add (g eta) to appropriate field
	 Vmath::Vadd(nq, flux[i+1][i], 1, tmp, 1, flux[i+1][i], 1);
       }

  }
  
  void LinearSWE::ConservativeToPrimitive(const Array<OneD, const Array<OneD, NekDouble> >&physin,
					      Array<OneD,       Array<OneD, NekDouble> >&physout)
  {
    int nq = GetTotPoints();
587
      
Claes Eskilsson's avatar
Claes Eskilsson committed
588
    if(physin.get() == physout.get())
589
      {
Claes Eskilsson's avatar
Claes Eskilsson committed
590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606
	// copy indata and work with tmp array
	Array<OneD, Array<OneD, NekDouble> >tmp(3);
	for (int i = 0; i < 3; ++i)
	  {
	    // deep copy
	    tmp[i] = Array<OneD, NekDouble>(nq);
	    Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
	  }
	
	// \eta = h - d
	Vmath::Vsub(nq,tmp[0],1,m_depth,1,physout[0],1);
	
	// u = hu/h
	Vmath::Vdiv(nq,tmp[1],1,tmp[0],1,physout[1],1);
	
	// v = hv/ v
	Vmath::Vdiv(nq,tmp[2],1,tmp[0],1,physout[2],1);
607
      }
Claes Eskilsson's avatar
Claes Eskilsson committed
608
    else
609
      {
Claes Eskilsson's avatar
Claes Eskilsson committed
610 611 612 613 614 615 616 617
	// \eta = h - d
	Vmath::Vsub(nq,physin[0],1,m_depth,1,physout[0],1);
	
	// u = hu/h
	Vmath::Vdiv(nq,physin[1],1,physin[0],1,physout[1],1);
	
	// v = hv/ v
	Vmath::Vdiv(nq,physin[2],1,physin[0],1,physout[2],1);
618 619
      }
  }
Claes Eskilsson's avatar
Claes Eskilsson committed
620 621 622


   void LinearSWE::v_ConservativeToPrimitive( )
623
  {
Claes Eskilsson's avatar
Claes Eskilsson committed
624 625 626 627 628 629 630
    int nq = GetTotPoints();
    
    // u = hu/h
    Vmath::Vdiv(nq,m_fields[1]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
	
    // v = hv/ v
    Vmath::Vdiv(nq,m_fields[2]->GetPhys(),1,m_fields[0]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
631

Claes Eskilsson's avatar
Claes Eskilsson committed
632 633
    // \eta = h - d
    Vmath::Vsub(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
634 635
  }

Claes Eskilsson's avatar
Claes Eskilsson committed
636 637
  void LinearSWE::PrimitiveToConservative(const Array<OneD, const Array<OneD, NekDouble> >&physin,
					     Array<OneD,       Array<OneD, NekDouble> >&physout)
638 639
  {
    
Claes Eskilsson's avatar
Claes Eskilsson committed
640
    int nq = GetTotPoints();
641
    
Claes Eskilsson's avatar
Claes Eskilsson committed
642
    if(physin.get() == physout.get())
643
      {
Claes Eskilsson's avatar
Claes Eskilsson committed
644 645 646
	// copy indata and work with tmp array
	Array<OneD, Array<OneD, NekDouble> >tmp(3);
	for (int i = 0; i < 3; ++i)
647
	  {
Claes Eskilsson's avatar
Claes Eskilsson committed
648 649 650
	    // deep copy
	    tmp[i] = Array<OneD, NekDouble>(nq);
	    Vmath::Vcopy(nq,physin[i],1,tmp[i],1);
651 652
	  }
	
Claes Eskilsson's avatar
Claes Eskilsson committed
653 654
	// h = \eta + d
	Vmath::Vadd(nq,tmp[0],1,m_depth,1,physout[0],1);
655
	
Claes Eskilsson's avatar
Claes Eskilsson committed
656 657
	// hu = h * u
	Vmath::Vmul(nq,physout[0],1,tmp[1],1,physout[1],1);
658
	
Claes Eskilsson's avatar
Claes Eskilsson committed
659 660 661
	// hv = h * v
	Vmath::Vmul(nq,physout[0],1,tmp[2],1,physout[2],1);
      
662
      }
Claes Eskilsson's avatar
Claes Eskilsson committed
663
    else
664
      {
Claes Eskilsson's avatar
Claes Eskilsson committed
665 666
	// h = \eta + d
	Vmath::Vadd(nq,physin[0],1,m_depth,1,physout[0],1);
667
	
Claes Eskilsson's avatar
Claes Eskilsson committed
668 669
	// hu = h * u
	Vmath::Vmul(nq,physout[0],1,physin[1],1,physout[1],1);
670
	
Claes Eskilsson's avatar
Claes Eskilsson committed
671 672 673 674 675
	// hv = h * v
	Vmath::Vmul(nq,physout[0],1,physin[2],1,physout[2],1);
	
      }
     
676 677
  }

Claes Eskilsson's avatar
Claes Eskilsson committed
678
  void LinearSWE::v_PrimitiveToConservative( )
679
  {
Claes Eskilsson's avatar
Claes Eskilsson committed
680
    int nq = GetTotPoints();
681
    
Claes Eskilsson's avatar
Claes Eskilsson committed
682 683
    // h = \eta + d
    Vmath::Vadd(nq,m_fields[0]->GetPhys(),1,m_depth,1,m_fields[0]->UpdatePhys(),1);
684
    
Claes Eskilsson's avatar
Claes Eskilsson committed
685 686
    // hu = h * u
    Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[1]->GetPhys(),1,m_fields[1]->UpdatePhys(),1);
687
    
Claes Eskilsson's avatar
Claes Eskilsson committed
688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712
    // hv = h * v
    Vmath::Vmul(nq,m_fields[0]->GetPhys(),1,m_fields[2]->GetPhys(),1,m_fields[2]->UpdatePhys(),1);
  }


 /**
     * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
     * \f$ h\mathbf{v} \f$.
     * 
     * @param physfield  Velocity field.
     * @param velocity   Velocity field.
     */
    void LinearSWE::GetVelocityVector(
        const Array<OneD, Array<OneD, NekDouble> > &physfield,
              Array<OneD, Array<OneD, NekDouble> > &velocity)
    {
        const int npts = physfield[0].num_elements();
        
        for (int i = 0; i < m_spacedim; ++i)
        {
            Vmath::Vcopy(npts, physfield[1+i], 1, velocity[i], 1);
        }
    }


713
    void LinearSWE::v_GenerateSummary(SolverUtils::SummaryList& s)
Claes Eskilsson's avatar
Claes Eskilsson committed
714
    {
715
        ShallowWaterSystem::v_GenerateSummary(s);
Claes Eskilsson's avatar
Claes Eskilsson committed
716 717 718 719 720 721
	if (m_session->DefinesSolverInfo("UpwindType"))
	  {
	    std::string UpwindType;
	    UpwindType = m_session->GetSolverInfo("UpwindType");
	    if (UpwindType == "LinearAverage")
	      {
722
	        SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear Average");
Claes Eskilsson's avatar
Claes Eskilsson committed
723 724 725
	      }
	    if (UpwindType == "LinearHLL")
	      {
726
	        SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear HLL");
Claes Eskilsson's avatar
Claes Eskilsson committed
727 728
	      }
	  }
729 730 731
	SolverUtils::AddSummaryItem(s, "Variables", "eta  should be in field[0]");
	SolverUtils::AddSummaryItem(s, "",          "u    should be in field[1]");
	SolverUtils::AddSummaryItem(s, "",          "v    should be in field[2]");
Claes Eskilsson's avatar
Claes Eskilsson committed
732 733
    }

734 735
} //end of namespace