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 using namespace std;
41 
42 namespace Nektar
43 {
44  string EulerADCFE::className =
46  "EulerADCFE", EulerADCFE::create,
47  "Euler equations in conservative variables with "
48  "artificial diffusion.");
49 
50  EulerADCFE::EulerADCFE(
52  : CompressibleFlowSystem(pSession)
53  {
54  }
55 
57  {
59 
60  if (m_shockCaptureType == "Smooth")
61  {
62  ASSERTL0(m_fields.num_elements() == m_spacedim + 3,
63  "Not enough variables for smooth shock capturing; "
64  "make sure you have added eps to variable list.");
65  m_smoothDiffusion = true;
66  }
67 
68  m_diffusion->SetArtificialDiffusionVector(
70 
71  if(m_session->DefinesSolverInfo("PROBLEMTYPE"))
72  {
73  std::string ProblemTypeStr = m_session->GetSolverInfo("PROBLEMTYPE");
74  int i;
75  for(i = 0; i < (int) SIZE_ProblemType; ++i)
76  {
77  if(boost::iequals(ProblemTypeMap[i], ProblemTypeStr))
78  {
80  break;
81  }
82  }
83  }
84  else
85  {
87  }
88 
90  {
92  DoOdeRhs, this);
94  DoOdeProjection, this);
95  }
96  else
97  {
98  ASSERTL0(false, "Implicit CFE not set up.");
99  }
100  }
101 
103  {
104 
105  }
106 
108  {
111  s, "Problem Type", ProblemTypeMap[m_problemType]);
112  }
113 
115  NekDouble initialtime,
116  bool dumpInitialConditions,
117  const int domain)
118  {
119  EquationSystem::v_SetInitialConditions(initialtime, false);
121 
122  if(dumpInitialConditions)
123  {
124  // Dump initial conditions to file
126  }
127  }
128 
130  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
131  Array<OneD, Array<OneD, NekDouble> > &outarray,
132  const NekDouble time)
133  {
134  int i;
135  int nvariables = inarray.num_elements();
136  int npoints = GetNpoints();
137 
139  Array<OneD, Array<OneD, NekDouble> > outarrayAdv(nvariables);
140  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
141 
142  for (i = 0; i < nvariables; ++i)
143  {
144  outarrayAdv[i] = Array<OneD, NekDouble>(npoints, 0.0);
145  outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
146  }
147 
148  m_advection->Advect(nvariables, m_fields, advVel, inarray,
149  outarrayAdv, m_time);
150 
151  for (i = 0; i < nvariables; ++i)
152  {
153  Vmath::Neg(npoints, outarrayAdv[i], 1);
154  }
155 
156  m_diffusion->Diffuse(nvariables, m_fields, inarray, outarrayDiff);
157 
158  if (m_shockCaptureType == "NonSmooth")
159  {
160  for (i = 0; i < nvariables; ++i)
161  {
162  Vmath::Vadd(npoints,
163  outarrayAdv[i], 1,
164  outarrayDiff[i], 1,
165  outarray[i], 1);
166  }
167  }
168  if(m_shockCaptureType == "Smooth")
169  {
170  const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();
171 
172  NekDouble pOrder = Vmath::Vmax(ExpOrder.num_elements(), ExpOrder, 1);
173 
174  Array <OneD, NekDouble > a_vel (npoints, 0.0);
175  Array <OneD, NekDouble > u_abs (npoints, 0.0);
176  Array <OneD, NekDouble > pres (npoints, 0.0);
177  Array <OneD, NekDouble > wave_sp(npoints, 0.0);
178 
179  GetPressure(inarray, pres);
180  GetSoundSpeed(inarray, pres, a_vel);
181  GetAbsoluteVelocity(inarray, u_abs);
182 
183  Vmath::Vadd(npoints, a_vel, 1, u_abs, 1, wave_sp, 1);
184 
185  NekDouble max_wave_sp = Vmath::Vmax(npoints, wave_sp, 1);
186 
187  Vmath::Smul(npoints,
188  m_C2,
189  outarrayDiff[nvariables-1], 1,
190  outarrayDiff[nvariables-1], 1);
191 
192  Vmath::Smul(npoints,
193  max_wave_sp,
194  outarrayDiff[nvariables-1], 1,
195  outarrayDiff[nvariables-1], 1);
196 
197  Vmath::Smul(npoints,
198  pOrder,
199  outarrayDiff[nvariables-1], 1,
200  outarrayDiff[nvariables-1], 1);
201 
202  for (i = 0; i < nvariables; ++i)
203  {
204  Vmath::Vadd(npoints,
205  outarrayAdv[i], 1,
206  outarrayDiff[i], 1,
207  outarray[i], 1);
208  }
209 
210  Array<OneD, Array<OneD, NekDouble> > outarrayForcing(nvariables);
211 
212  for (i = 0; i < nvariables; ++i)
213  {
214  outarrayForcing[i] = Array<OneD, NekDouble>(npoints, 0.0);
215  }
216 
217  GetForcingTerm(inarray, outarrayForcing);
218 
219  for (i = 0; i < nvariables; ++i)
220  {
221  // Add Forcing Term
222  Vmath::Vadd(npoints,
223  outarray[i], 1,
224  outarrayForcing[i], 1,
225  outarray[i], 1);
226  }
227  }
228 
229  // Add sponge layer if defined in the session file
230  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
231  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
232  {
233  (*x)->Apply(m_fields, inarray, outarray, time);
234  }
235  }
236 
238  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
239  Array<OneD, Array<OneD, NekDouble> > &outarray,
240  const NekDouble time)
241  {
242  int i;
243  int nvariables = inarray.num_elements();
244 
245  switch(m_projectionType)
246  {
248  {
249  // Just copy over array
250  int npoints = GetNpoints();
251 
252  for(i = 0; i < nvariables; ++i)
253  {
254  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
255  }
256  SetBoundaryConditions(outarray, time);
257  break;
258  }
261  {
262  ASSERTL0(false, "No Continuous Galerkin for Euler equations");
263  break;
264  }
265  default:
266  ASSERTL0(false, "Unknown projection scheme");
267  break;
268  }
269  }
270 
272  Array<OneD, Array<OneD, NekDouble> > &inarray,
273  NekDouble time)
274  {
275  std::string varName;
276  int cnt = 0;
277 
278  // loop over Boundary Regions
279  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
280  {
281  std::string type = m_fields[0]->GetBndConditions()[n]->GetUserDefined();
282 
283  // Wall Boundary Condition
284  if (boost::iequals(type,"WallViscous"))
285  {
286  // Wall Boundary Condition
287  ASSERTL0(false, "WallViscous is a wrong bc for the "
288  "Euler equations");
289  }
290  else
291  {
292  SetCommonBC(type,n,time, cnt,inarray);
293  }
294 
295  // no User Defined conditions provided so skip cnt
296  // this line is left in case solver specific condition is added.
297  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
298  }
299  }
300 }
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Definition: EulerADCFE.cpp:271
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Definition: EulerADCFE.cpp:114
void GetSoundSpeed(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &soundspeed)
NekDouble m_time
Current time of simulation.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp()
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
virtual ~EulerADCFE()
problem type selector
Definition: EulerADCFE.cpp:102
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
STL namespace.
void GetArtificialDynamicViscosity(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &mu_var)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
ProblemType m_problemType
Definition: EulerADCFE.h:79
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition: EulerADCFE.cpp:129
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition: EulerADCFE.cpp:237
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
Definition: EulerADCFE.cpp:56
void SetCommonBC(const std::string &userDefStr, const int n, const NekDouble time, int &cnt, Array< OneD, Array< OneD, NekDouble > > &inarray)
Set boundary conditions which can be: a) Wall and Symmerty BCs implemented at CompressibleFlowSystem ...
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
Definition: EulerADCFE.cpp:107
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
int m_spacedim
Spatial dimension (>= expansion dim).
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Length of enum list.
Definition: EulerADCFE.h:47
EquationSystemFactory & GetEquationSystemFactory()
void GetAbsoluteVelocity(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &Vtot)
SolverUtils::AdvectionSharedPtr m_advection
SolverUtils::DiffusionSharedPtr m_diffusion
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
void GetForcingTerm(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > outarrayForcing)
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
ProblemType
Definition: EulerADCFE.h:44
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
const char *const ProblemTypeMap[]
Definition: EulerADCFE.h:50