Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
APEUpwindSolver.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: APEUpwindSolver.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: Upwind Riemann solver for the APE equations.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42 
43 std::string APEUpwindSolver::solverName = SolverUtils::GetRiemannSolverFactory().
44  RegisterCreatorFunction("APEUpwind", APEUpwindSolver::create, "APE Upwind solver");
45 
46 /**
47 *
48 */
49 
50 /**
51 *
52 */
53 APEUpwindSolver::APEUpwindSolver() :
55 {
56  // we need out own rotation logic
57  m_requiresRotation = false;
58 }
59 
60 /**
61 *
62 */
64  const int nDim,
65  const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
66  const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
67  Array<OneD, Array<OneD, NekDouble> > &flux)
68 {
69 
70  ASSERTL1(CheckVectors("N"), "N not defined.");
71  ASSERTL1(CheckAuxVec ("vecLocs"), "vecLoc not defined.");
72  ASSERTL1(CheckVectors("basefield"), "basefield not defined.");
73  const Array<OneD, const Array<OneD, NekDouble> > normals =
74  m_vectors["N"]();
75  const Array<OneD, const Array<OneD, NekDouble> > vecLocs =
76  m_auxVec["vecLocs"]();
77  const Array<OneD, const Array<OneD, NekDouble> > basefield =
78  m_vectors["basefield"]();
79 
80  int nFields = Fwd .num_elements();
81  int nPts = Fwd[0].num_elements();
82 
83  // rotate and store basefield
84  m_rotBasefield = Array<OneD, Array<OneD, NekDouble> > (nDim+1);
85  for (int i = 0; i < nDim + 1; i++)
86  {
87  m_rotBasefield[i] = Array<OneD, NekDouble>(nPts);
88  }
89  Array<OneD, Array<OneD, NekDouble> > baseVecLocs(1);
90  baseVecLocs[0] = Array<OneD, NekDouble>(nDim);
91  for (int i = 0; i < nDim; ++i)
92  {
93  baseVecLocs[0][i] = 1+i;
94  }
95 
96  rotateToNormal(basefield, normals, baseVecLocs, m_rotBasefield);
97 
98  if (m_rotStorage[0].num_elements() != nFields ||
99  m_rotStorage[0][0].num_elements() != nPts)
100  {
101  for (int i = 0; i < 3; ++i)
102  {
103  m_rotStorage[i] =
104  Array<OneD, Array<OneD, NekDouble> >(nFields);
105  for (int j = 0; j < nFields; ++j)
106  {
107  m_rotStorage[i][j] = Array<OneD, NekDouble>(nPts);
108  }
109  }
110  }
111 
112  rotateToNormal(Fwd, normals, vecLocs, m_rotStorage[0]);
113  rotateToNormal(Bwd, normals, vecLocs, m_rotStorage[1]);
115  rotateFromNormal(m_rotStorage[2], normals, vecLocs, flux);
116 }
117 
118 
119 
120 /**
121 *
122 */
124  const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
125  const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
126  Array<OneD, Array<OneD, NekDouble> > &flux)
127 {
128  // fetch params
129  ASSERTL1(CheckParams("Gamma"), "Gamma not defined.");
130  ASSERTL1(CheckParams("Rho"), "Rho not defined.");
131  ASSERTL1(CheckVectors("N"), "N not defined.");
132  const NekDouble &gamma= m_params["Gamma"]();
133  const NekDouble &rho = m_params["Rho"]();
134  const Array<OneD, const Array<OneD, NekDouble> > normals = m_vectors["N"]();
135  const Array<OneD, const Array<OneD, NekDouble> > basefield = m_rotBasefield;
136 
137  int dim = normals.num_elements();
138  int nvar = dim +1;
139  ASSERTL1(Fwd.num_elements() == nvar, "Fwd malformed.");
140  ASSERTL1(Bwd.num_elements() == nvar, "Bwd malformed.");
141 
142  int nTracePts = Fwd[0].num_elements();
143 
144  Array<OneD, Array<OneD, NekDouble> > upphysfield(2);
145  for (int i = 0; i < 2; i++)
146  {
147  upphysfield[i] = Array<OneD, NekDouble>(nTracePts);
148  }
149 
150  for (int i = 0; i < nTracePts; i++)
151  {
152  // assign variables
153  NekDouble pL = Fwd[0][i];
154  NekDouble uL = Fwd[1][i];
155 
156  NekDouble pR = Bwd[0][i];
157  NekDouble uR = Bwd[1][i];
158 
159  NekDouble p0 = basefield[0][i];
160  NekDouble u0 = basefield[1][i];
161 
162  Array<OneD, NekDouble> characteristic(4);
163  Array<OneD, NekDouble> W(2);
164  Array<OneD, NekDouble> lambda(2);
165 
166  // compute the wave speeds
167  lambda[0]=u0 + sqrt(p0*gamma*rho)/rho;
168  lambda[1]=u0 - sqrt(p0*gamma*rho)/rho;
169 
170  // calculate the caracteristic variables
171  //left characteristics
172  characteristic[0] = pL/2 + uL*sqrt(p0*gamma*rho)/2;
173  characteristic[1] = pL/2 - uL*sqrt(p0*gamma*rho)/2;
174  //right characteristics
175  characteristic[2] = pR/2 + uR*sqrt(p0*gamma*rho)/2;
176  characteristic[3] = pR/2 - uR*sqrt(p0*gamma*rho)/2;
177 
178  //take left or right value of characteristic variable
179  for (int j = 0; j < 2; j++)
180  {
181  if (lambda[j]>=0)
182  {
183  W[j]=characteristic[j];
184  }
185  if(lambda[j]<0)
186  {
187  W[j]=characteristic[j+2];
188  }
189  }
190 
191  //calculate conservative variables from characteristics
192  upphysfield[0][i] = W[0]+W[1];
193  // do not divide by zero
194  if (p0*gamma*rho != 0)
195  {
196  upphysfield[1][i] = (W[0]-W[1])/sqrt(p0*gamma*rho);
197  }
198  else
199  {
200  upphysfield[1][i] = 0.0;
201  }
202  }
203 
204  // compute the fluxes
205 
206  // flux[0][i] = gamma*p0 * upphysfield[1] + u0*upphysfield[0]
207  Vmath::Smul(nTracePts, gamma, basefield[0], 1, flux[0], 1);
208  Vmath::Vmul(nTracePts, upphysfield[1], 1, flux[0], 1, flux[0], 1);
209  Vmath::Vvtvp(nTracePts, basefield[1], 1, upphysfield[0], 1, flux[0], 1, flux[0], 1);
210 
211  // flux[1][i] = upphysfield[0]/rho + u0*upphysfield[1];
212  Vmath::Smul(nTracePts, 1/rho, upphysfield[0], 1, flux[1], 1);
213  Vmath::Vvtvp(nTracePts, basefield[1], 1, upphysfield[1], 1, flux[1], 1, flux[1], 1);
214  for (int j = 2; j < nvar; j++)
215  {
216  // flux[1][i] = upphysfield[0]/rho + u0*upphysfield[1] + v0 * vL + w0 + wL;
217  Vmath::Vvtvp(nTracePts, basefield[j], 1, Fwd[j], 1, flux[1], 1, flux[1], 1);
218 
219  // flux[2/3][i] = 0.0;
220  Vmath::Zero(nTracePts, flux[j], 1);
221  }
222 }
223 }