Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
EulerADCFE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File EulerArtificialDiffusionCFE.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: Euler equations in conservative variables with artificial
33 // diffusion
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
38 #include <boost/algorithm/string.hpp>
39 
40 namespace Nektar
41 {
42  string EulerADCFE::className =
44  "EulerADCFE", EulerADCFE::create,
45  "Euler equations in conservative variables with "
46  "artificial diffusion.");
47 
50  : CompressibleFlowSystem(pSession)
51  {
52  }
53 
55  {
57 
58  if (m_shockCaptureType == "Smooth")
59  {
60  ASSERTL0(m_fields.num_elements() == m_spacedim + 3,
61  "Not enough variables for smooth shock capturing; "
62  "make sure you have added eps to variable list.");
63  m_smoothDiffusion = true;
64  }
65 
66  m_diffusion->SetArtificialDiffusionVector(
68 
69  if(m_session->DefinesSolverInfo("PROBLEMTYPE"))
70  {
71  std::string ProblemTypeStr = m_session->GetSolverInfo("PROBLEMTYPE");
72  int i;
73  for(i = 0; i < (int) SIZE_ProblemType; ++i)
74  {
75  if(boost::iequals(ProblemTypeMap[i], ProblemTypeStr))
76  {
78  break;
79  }
80  }
81  }
82  else
83  {
85  }
86 
88  {
90  DoOdeRhs, this);
92  DoOdeProjection, this);
93  }
94  else
95  {
96  ASSERTL0(false, "Implicit CFE not set up.");
97  }
98  }
99 
101  {
102 
103  }
104 
106  {
109  s, "Problem Type", ProblemTypeMap[m_problemType]);
110  }
111 
113  NekDouble initialtime,
114  bool dumpInitialConditions,
115  const int domain)
116  {
117  EquationSystem::v_SetInitialConditions(initialtime, false);
118 
119  if(dumpInitialConditions)
120  {
121  // Dump initial conditions to file
123  }
124  }
125 
127  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
128  Array<OneD, Array<OneD, NekDouble> > &outarray,
129  const NekDouble time)
130  {
131  int i;
132  int nvariables = inarray.num_elements();
133  int npoints = GetNpoints();
134 
135  Array<OneD, Array<OneD, NekDouble> > advVel;
136  Array<OneD, Array<OneD, NekDouble> > outarrayAdv(nvariables);
137  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
138 
139  for (i = 0; i < nvariables; ++i)
140  {
141  outarrayAdv[i] = Array<OneD, NekDouble>(npoints, 0.0);
142  outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
143  }
144 
145  m_advection->Advect(nvariables, m_fields, advVel, inarray, outarrayAdv);
146 
147  for (i = 0; i < nvariables; ++i)
148  {
149  Vmath::Neg(npoints, outarrayAdv[i], 1);
150  }
151 
152  m_diffusion->Diffuse(nvariables, m_fields, inarray, outarrayDiff);
153 
154  if (m_shockCaptureType == "NonSmooth")
155  {
156  for (i = 0; i < nvariables; ++i)
157  {
158  Vmath::Vadd(npoints,
159  outarrayAdv[i], 1,
160  outarrayDiff[i], 1,
161  outarray[i], 1);
162  }
163  }
164  if(m_shockCaptureType == "Smooth")
165  {
166  const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();
167 
168  NekDouble pOrder = Vmath::Vmax(ExpOrder.num_elements(), ExpOrder, 1);
169 
170  Array <OneD, NekDouble > a_vel (npoints, 0.0);
171  Array <OneD, NekDouble > u_abs (npoints, 0.0);
172  Array <OneD, NekDouble > pres (npoints, 0.0);
173  Array <OneD, NekDouble > wave_sp(npoints, 0.0);
174 
175  GetPressure(inarray, pres);
176  GetSoundSpeed(inarray, pres, a_vel);
177  GetAbsoluteVelocity(inarray, u_abs);
178 
179  Vmath::Vadd(npoints, a_vel, 1, u_abs, 1, wave_sp, 1);
180 
181  NekDouble max_wave_sp = Vmath::Vmax(npoints, wave_sp, 1);
182 
183  Vmath::Smul(npoints,
184  m_C2,
185  outarrayDiff[nvariables-1], 1,
186  outarrayDiff[nvariables-1], 1);
187 
188  Vmath::Smul(npoints,
189  max_wave_sp,
190  outarrayDiff[nvariables-1], 1,
191  outarrayDiff[nvariables-1], 1);
192 
193  Vmath::Smul(npoints,
194  pOrder,
195  outarrayDiff[nvariables-1], 1,
196  outarrayDiff[nvariables-1], 1);
197 
198  for (i = 0; i < nvariables; ++i)
199  {
200  Vmath::Vadd(npoints,
201  outarrayAdv[i], 1,
202  outarrayDiff[i], 1,
203  outarray[i], 1);
204  }
205 
206  Array<OneD, Array<OneD, NekDouble> > outarrayForcing(nvariables);
207 
208  for (i = 0; i < nvariables; ++i)
209  {
210  outarrayForcing[i] = Array<OneD, NekDouble>(npoints, 0.0);
211  }
212 
213  GetForcingTerm(inarray, outarrayForcing);
214 
215  for (i = 0; i < nvariables; ++i)
216  {
217  // Add Forcing Term
218  Vmath::Vadd(npoints,
219  outarray[i], 1,
220  outarrayForcing[i], 1,
221  outarray[i], 1);
222  }
223  }
224  }
225 
227  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
228  Array<OneD, Array<OneD, NekDouble> > &outarray,
229  const NekDouble time)
230  {
231  int i;
232  int nvariables = inarray.num_elements();
233 
234  switch(m_projectionType)
235  {
237  {
238  // Just copy over array
239  int npoints = GetNpoints();
240 
241  for(i = 0; i < nvariables; ++i)
242  {
243  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
244  }
245  SetBoundaryConditions(outarray, time);
246  break;
247  }
250  {
251  ASSERTL0(false, "No Continuous Galerkin for Euler equations");
252  break;
253  }
254  default:
255  ASSERTL0(false, "Unknown projection scheme");
256  break;
257  }
258  }
259 
261  Array<OneD, Array<OneD, NekDouble> > &inarray,
262  NekDouble time)
263  {
264  std::string varName;
265  int nvariables = m_fields.num_elements();
266  int cnt = 0;
267 
268  // loop over Boundary Regions
269  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
270  {
271  // Wall Boundary Condition
272  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
274  {
275  WallBC(n, cnt, inarray);
276  }
277 
278  // Wall Boundary Condition
279  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
281  {
282  ASSERTL0(false, "WallViscous is a wrong bc for the "
283  "Euler equations");
284  }
285 
286  // Symmetric Boundary Condition
287  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
289  {
290  SymmetryBC(n, cnt, inarray);
291  }
292 
293  // Riemann invariant characteristic Boundary Condition (CBC)
294  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
296  {
297  RiemannInvariantBC(n, cnt, inarray);
298  }
299 
300  // Extrapolation of the data at the boundaries
301  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
303  {
304  ExtrapOrder0BC(n, cnt, inarray);
305  }
306 
307  // Time Dependent Boundary Condition (specified in meshfile)
308  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined()
310  {
311  for (int i = 0; i < nvariables; ++i)
312  {
313  varName = m_session->GetVariable(i);
314  m_fields[i]->EvaluateBoundaryConditions(time, varName);
315  }
316  }
317 
318  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
319  }
320  }
321 }