Nektar++
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);
119 
120  if(dumpInitialConditions)
121  {
122  // Dump initial conditions to file
124  }
125  }
126 
128  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
129  Array<OneD, Array<OneD, NekDouble> > &outarray,
130  const NekDouble time)
131  {
132  int i;
133  int nvariables = inarray.num_elements();
134  int npoints = GetNpoints();
135 
137  Array<OneD, Array<OneD, NekDouble> > outarrayAdv(nvariables);
138  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
139 
140  for (i = 0; i < nvariables; ++i)
141  {
142  outarrayAdv[i] = Array<OneD, NekDouble>(npoints, 0.0);
143  outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
144  }
145 
146  m_advection->Advect(nvariables, m_fields, advVel, inarray,
147  outarrayAdv, m_time);
148 
149  for (i = 0; i < nvariables; ++i)
150  {
151  Vmath::Neg(npoints, outarrayAdv[i], 1);
152  }
153 
154  m_diffusion->Diffuse(nvariables, m_fields, inarray, outarrayDiff);
155 
156  if (m_shockCaptureType == "NonSmooth")
157  {
158  for (i = 0; i < nvariables; ++i)
159  {
160  Vmath::Vadd(npoints,
161  outarrayAdv[i], 1,
162  outarrayDiff[i], 1,
163  outarray[i], 1);
164  }
165  }
166  if(m_shockCaptureType == "Smooth")
167  {
168  const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();
169 
170  NekDouble pOrder = Vmath::Vmax(ExpOrder.num_elements(), ExpOrder, 1);
171 
172  Array <OneD, NekDouble > a_vel (npoints, 0.0);
173  Array <OneD, NekDouble > u_abs (npoints, 0.0);
174  Array <OneD, NekDouble > pres (npoints, 0.0);
175  Array <OneD, NekDouble > wave_sp(npoints, 0.0);
176 
177  GetPressure(inarray, pres);
178  GetSoundSpeed(inarray, pres, a_vel);
179  GetAbsoluteVelocity(inarray, u_abs);
180 
181  Vmath::Vadd(npoints, a_vel, 1, u_abs, 1, wave_sp, 1);
182 
183  NekDouble max_wave_sp = Vmath::Vmax(npoints, wave_sp, 1);
184 
185  Vmath::Smul(npoints,
186  m_C2,
187  outarrayDiff[nvariables-1], 1,
188  outarrayDiff[nvariables-1], 1);
189 
190  Vmath::Smul(npoints,
191  max_wave_sp,
192  outarrayDiff[nvariables-1], 1,
193  outarrayDiff[nvariables-1], 1);
194 
195  Vmath::Smul(npoints,
196  pOrder,
197  outarrayDiff[nvariables-1], 1,
198  outarrayDiff[nvariables-1], 1);
199 
200  for (i = 0; i < nvariables; ++i)
201  {
202  Vmath::Vadd(npoints,
203  outarrayAdv[i], 1,
204  outarrayDiff[i], 1,
205  outarray[i], 1);
206  }
207 
208  Array<OneD, Array<OneD, NekDouble> > outarrayForcing(nvariables);
209 
210  for (i = 0; i < nvariables; ++i)
211  {
212  outarrayForcing[i] = Array<OneD, NekDouble>(npoints, 0.0);
213  }
214 
215  GetForcingTerm(inarray, outarrayForcing);
216 
217  for (i = 0; i < nvariables; ++i)
218  {
219  // Add Forcing Term
220  Vmath::Vadd(npoints,
221  outarray[i], 1,
222  outarrayForcing[i], 1,
223  outarray[i], 1);
224  }
225  }
226 
227  // Add sponge layer if defined in the session file
228  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
229  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
230  {
231  (*x)->Apply(m_fields, inarray, outarray, time);
232  }
233  }
234 
236  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
237  Array<OneD, Array<OneD, NekDouble> > &outarray,
238  const NekDouble time)
239  {
240  int i;
241  int nvariables = inarray.num_elements();
242 
243  switch(m_projectionType)
244  {
246  {
247  // Just copy over array
248  int npoints = GetNpoints();
249 
250  for(i = 0; i < nvariables; ++i)
251  {
252  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
253  }
254  SetBoundaryConditions(outarray, time);
255  break;
256  }
259  {
260  ASSERTL0(false, "No Continuous Galerkin for Euler equations");
261  break;
262  }
263  default:
264  ASSERTL0(false, "Unknown projection scheme");
265  break;
266  }
267  }
268 
270  Array<OneD, Array<OneD, NekDouble> > &inarray,
271  NekDouble time)
272  {
273  std::string varName;
274  int cnt = 0;
275 
276  // loop over Boundary Regions
277  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
278  {
279  std::string type = m_fields[0]->GetBndConditions()[n]->GetUserDefined();
280 
281  // Wall Boundary Condition
282  if (boost::iequals(type,"WallViscous"))
283  {
284  // Wall Boundary Condition
285  ASSERTL0(false, "WallViscous is a wrong bc for the "
286  "Euler equations");
287  }
288  else
289  {
290  SetCommonBC(type,n,time, cnt,inarray);
291  }
292 
293  // no User Defined conditions provided so skip cnt
294  // this line is left in case solver specific condition is added.
295  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
296  }
297  }
298 }
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Definition: EulerADCFE.cpp:269
#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:112
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:756
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:100
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
void GetArtificialDynamicViscosity(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &mu_var)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
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:127
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:235
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
Definition: EulerADCFE.cpp:54
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:105
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
Length of enum list.
Definition: EulerADCFE.h:47
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
Definition: EulerADCFE.h:65
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)
static std::string className
Name of class.
Definition: EulerADCFE.h:74
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
EulerADCFE(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: EulerADCFE.cpp:48
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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