Nektar++
CFLtester.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File UnsteadyAdvection.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: CFL tester solve routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 
39 #include <LocalRegions/TriExp.h>
40 #include <LocalRegions/QuadExp.h>
43 
44 namespace Nektar
45 {
46  string CFLtester::className = GetEquationSystemFactory().RegisterCreatorFunction("CFLtester", CFLtester::create, "Testing CFL restriction");
47 
50  : UnsteadySystem(pSession)
51  {
52  }
53 
55  {
56  UnsteadySystem::v_InitObject();
57 
59  std::vector<std::string> vel;
60  vel.push_back("Vx");
61  vel.push_back("Vy");
62  vel.push_back("Vz");
63  vel.resize(m_spacedim);
64 
65  EvaluateFunction(vel, m_velocity, "AdvectionVelocity");
66 
68  {
71  }
72  else
73  {
74  ASSERTL0(false, "Implicit unsteady Advection not set up.");
75  }
76  }
77 
79  {
80  }
81 
83  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
84  Array<OneD, Array<OneD, NekDouble> >&outarray,
85  const NekDouble time)
86  {
87  int j;
88  int nvariables = inarray.num_elements();
89  int npoints = GetNpoints();
90 
91  switch (m_projectionType)
92  {
94  {
95  int ncoeffs = inarray[0].num_elements();
96  Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
97 
98  WeakAdv[0] = Array<OneD, NekDouble>(ncoeffs*nvariables);
99  for(j = 1; j < nvariables; ++j)
100  {
101  WeakAdv[j] = WeakAdv[j-1] + ncoeffs;
102  }
103 
104  WeakDGAdvection(inarray, WeakAdv, true, true);
105 
106  for(j = 0; j < nvariables; ++j)
107  {
108  m_fields[j]->MultiplyByElmtInvMass(WeakAdv[j], WeakAdv[j]);
109  m_fields[j]->BwdTrans(WeakAdv[j],outarray[j]);
110  Vmath::Neg(npoints,outarray[j],1);
111  }
112  break;
113  }
116  {
117  // Calculate -V\cdot Grad(u);
118  for(j = 0; j < nvariables; ++j)
119  {
121  inarray[j],
122  outarray[j]);
123 
124  Vmath::Neg(npoints, outarray[j], 1);
125  }
126  break;
127  }
128  }
129  }
130 
132  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
133  Array<OneD, Array<OneD, NekDouble> >&outarray,
134  const NekDouble time)
135  {
136  int j;
137  int nvariables = inarray.num_elements();
138  SetBoundaryConditions(time);
139 
140  switch(m_projectionType)
141  {
143  {
144  // Just copy over array
145  int npoints = GetNpoints();
146 
147  for(j = 0; j < nvariables; ++j)
148  {
149  Vmath::Vcopy(npoints,inarray[j],1,outarray[j],1);
150  }
151  }
152  break;
155  {
157 
158  for(j = 0; j < nvariables; ++j)
159  {
160  m_fields[j]->FwdTrans(inarray[j],coeffs);
161  m_fields[j]->BwdTrans_IterPerExp(coeffs,outarray[j]);
162  }
163  break;
164  }
165  default:
166  ASSERTL0(false,"Unknown projection scheme");
167  break;
168  }
169  }
170 
171 
173  const int i,
174  Array<OneD, Array<OneD, NekDouble> > &physfield,
176  {
177  ASSERTL1(flux.num_elements() == m_velocity.num_elements(),
178  "Dimension of flux array and velocity array do not match");
179 
180  for (int j = 0; j < flux.num_elements(); ++j)
181  {
183  physfield[i], 1,
184  m_velocity[j], 1,
185  flux[j], 1);
186  }
187  }
188 
190  Array<OneD, Array<OneD, NekDouble> > &physfield,
191  Array<OneD, Array<OneD, NekDouble> > &numflux)
192  {
193  int i;
194  int nTraceNumPoints = GetTraceNpoints();
195  int nvel = m_spacedim;
196 
197  Array<OneD, NekDouble > Fwd(nTraceNumPoints);
198  Array<OneD, NekDouble > Bwd(nTraceNumPoints);
199  Array<OneD, NekDouble > Vn (nTraceNumPoints,0.0);
200 
201  // Get Edge Velocity - Could be stored if time independent
202  for (i = 0; i < nvel; ++i)
203  {
204  m_fields[0]->ExtractTracePhys(m_velocity[i], Fwd);
205  Vmath::Vvtvp(nTraceNumPoints,
206  m_traceNormals[i], 1,
207  Fwd, 1, Vn, 1, Vn, 1);
208  }
209 
210  for (i = 0; i < numflux.num_elements(); ++i)
211  {
212  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
213  //evaulate upwinded m_fields[i]
214  m_fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
215  // calculate m_fields[i]*Vn
216  Vmath::Vmul(nTraceNumPoints, numflux[i], 1, Vn, 1, numflux[i], 1);
217  }
218  }
219 
220 
221 
223  {
224  UnsteadySystem::v_GenerateSummary(s);
225  }
226 
227 
228 
230  const Array<OneD,int> ExpOrder,
231  const Array<OneD,NekDouble> CFL,
232  NekDouble timeCFL)
233  {
234 
235  int n_element = m_fields[0]->GetExpSize();
236 
237  //const NekDouble minLengthStdTri = 1.414213;
238  //const NekDouble minLengthStdQuad = 2.0;
239  //const NekDouble cLambda = 0.2; // Spencer book pag. 317
240 
241  Array<OneD, NekDouble> tstep (n_element, 0.0);
242  Array<OneD, NekDouble> stdVelocity(n_element, 0.0);
243  stdVelocity = GetStdVelocity(m_velocity);
244 
245  for (int el = 0; el < n_element; ++el)
246  {
247  int npoints = m_fields[0]->GetExp(el)->GetTotPoints();
248  Array<OneD, NekDouble> one2D(npoints, 1.0);
249  //NekDouble Area = m_fields[0]->GetExp(el)->Integral(one2D);
250  if(boost::dynamic_pointer_cast<LocalRegions::TriExp>(m_fields[0]->GetExp(el)))
251  {
252  //tstep[el] = timeCFL/(stdVelocity[el]*cLambda*(ExpOrder[el]-1)*(ExpOrder[el]-1));
253  //tstep[el] = timeCFL*minLengthStdTri/(stdVelocity[el]*cLambda*(ExpOrder[el]-1)*(ExpOrder[el]-1));
254  tstep[el] = CFL[el]/(stdVelocity[el]);
255  }
256  else if(boost::dynamic_pointer_cast<LocalRegions::QuadExp>(m_fields[0]->GetExp(el)))
257  {
258  //tstep[el] = timeCFL/(stdVelocity[el]*cLambda*(ExpOrder[el]-1)*(ExpOrder[el]-1));
259  //tstep[el] = timeCFL*minLengthStdQuad/(stdVelocity[el]*cLambda*(ExpOrder[el]-1)*(ExpOrder[el]-1));
260  tstep[el] = CFL[el]/(stdVelocity[el]);
261  }
262  }
263 
264  NekDouble TimeStep = Vmath::Vmin(n_element, tstep, 1);
265 
266  return TimeStep;
267  }
268 
269 
270 
272  int ExpOrder,
273  NekDouble CFL,
274  NekDouble TimeStability)
275  {
276  //================================================================
277  // This function has been created just to test specific problems, hence is not general
278  // and it has been implemented in a rude fashion, as the full CFLtester class.
279  // For real CFL calculations refer to the general implementation above. (A.Bolis)
280  //================================================================
281 
282  NekDouble TimeStep;
283  NekDouble SpatialStability;
284  int n_elements = m_fields[0]->GetExpSize();
285 
286  //solve ambiguity in windows
287  NekDouble n_elem = n_elements;
288  NekDouble DH = sqrt(n_elem);
289 
290  int H = (int)DH;
291  int P = ExpOrder-1;
292 
293  //================================================================
294  // Regular meshes
295 
296  //SpatialStability = EigenvaluesRegMeshes[H-1][P-1];
297 
298  //================================================================
299  // Anisotropic meshes
300 
301  if (TimeStability == 1.0)
302  {
303  SpatialStability = EigenvaluesAnaMeshesAB2[H/2][P-1];
304  }
305  else if (TimeStability == 2.0)
306  {
307  SpatialStability = EigenvaluesAnaMeshesRK2[H/2][P-1];
308  }
309  else if (TimeStability == 2.784)
310  {
311  SpatialStability = EigenvaluesAnaMeshesRK4[H/2][P-1];
312  }
313  else
314  {
315  ASSERTL0(false,"Dominant eigenvalues database not present for this time-stepping scheme")
316  }
317 
318  //================================================================
319 
320  TimeStep = (TimeStability/SpatialStability)*CFL;
321 
322  //================================================================
323 
324  return TimeStep;
325  }
326 
327 
328 
330  const Array<OneD, Array<OneD,NekDouble> > inarray)
331  {
332  // Checking if the problem is 2D
333  ASSERTL0(m_expdim >= 2, "Method not implemented for 1D");
334 
335  int nTotQuadPoints = GetTotPoints();
336  int n_element = m_fields[0]->GetExpSize();
337  int nvel = inarray.num_elements();
338 
339  NekDouble pntVelocity;
340 
341  // Getting the standard velocity vector on the 2D normal space
342  Array<OneD, Array<OneD, NekDouble> > stdVelocity(nvel);
343  Array<OneD, NekDouble> stdV(n_element, 0.0);
344 
345  for (int i = 0; i < nvel; ++i)
346  {
347  stdVelocity[i] = Array<OneD, NekDouble>(nTotQuadPoints);
348  }
349 
350  if (nvel == 2)
351  {
352  for (int el = 0; el < n_element; ++el)
353  {
355  m_fields[0]->GetExp(el)->as<LocalRegions::Expansion2D>();
356  int n_points = el2D->GetTotPoints();
357 
359  el2D->GetGeom2D()->GetMetricInfo()->GetJac();
361  el2D->GetGeom2D()->GetMetricInfo()->GetDerivFactors();
362 
363  if (el2D->GetGeom2D()->GetMetricInfo()->GetGtype()
365  {
366  for (int i = 0; i < n_points; i++)
367  {
368  stdVelocity[0][i] = gmat[0][i]*inarray[0][i]
369  + gmat[2][i]*inarray[1][i];
370 
371  stdVelocity[1][i] = gmat[1][i]*inarray[0][i]
372  + gmat[3][i]*inarray[1][i];
373  }
374  }
375  else
376  {
377  for (int i = 0; i < n_points; i++)
378  {
379  stdVelocity[0][i] = gmat[0][0]*inarray[0][i]
380  + gmat[2][0]*inarray[1][i];
381 
382  stdVelocity[1][i] = gmat[1][0]*inarray[0][i]
383  + gmat[3][0]*inarray[1][i];
384  }
385  }
386 
387 
388  for (int i = 0; i < n_points; i++)
389  {
390  pntVelocity = sqrt(stdVelocity[0][i]*stdVelocity[0][i]
391  + stdVelocity[1][i]*stdVelocity[1][i]);
392 
393  if (pntVelocity>stdV[el])
394  {
395  stdV[el] = pntVelocity;
396  }
397  }
398  }
399  }
400  else
401  {
402  for (int el = 0; el < n_element; ++el)
403  {
405  m_fields[0]->GetExp(el)->as<LocalRegions::Expansion3D>();
406 
407  int n_points = el3D->GetTotPoints();
408 
410  el3D->GetGeom3D()->GetMetricInfo()->GetJac();
412  el3D->GetGeom3D()->GetMetricInfo()->GetDerivFactors();
413 
414  if (el3D->GetGeom3D()->GetMetricInfo()->GetGtype()
416  {
417  for (int i = 0; i < n_points; i++)
418  {
419  stdVelocity[0][i] = gmat[0][i]*inarray[0][i]
420  + gmat[3][i]*inarray[1][i]
421  + gmat[6][i]*inarray[2][i];
422 
423  stdVelocity[1][i] = gmat[1][i]*inarray[0][i]
424  + gmat[4][i]*inarray[1][i]
425  + gmat[7][i]*inarray[2][i];
426 
427  stdVelocity[2][i] = gmat[2][i]*inarray[0][i]
428  + gmat[5][i]*inarray[1][i]
429  + gmat[8][i]*inarray[2][i];
430  }
431  }
432  else
433  {
435  el3D->GetGeom3D()->GetMetricInfo()->GetJac();
437  el3D->GetGeom3D()->GetMetricInfo()->GetDerivFactors();
438 
439  for (int i = 0; i < n_points; i++)
440  {
441  stdVelocity[0][i] = gmat[0][0]*inarray[0][i]
442  + gmat[3][0]*inarray[1][i]
443  + gmat[6][0]*inarray[2][i];
444 
445  stdVelocity[1][i] = gmat[1][0]*inarray[0][i]
446  + gmat[4][0]*inarray[1][i]
447  + gmat[7][0]*inarray[2][i];
448 
449  stdVelocity[2][i] = gmat[2][0]*inarray[0][i]
450  + gmat[5][0]*inarray[1][i]
451  + gmat[8][0]*inarray[2][i];
452  }
453  }
454 
455  for (int i = 0; i < n_points; i++)
456  {
457  pntVelocity = sqrt(stdVelocity[0][i]*stdVelocity[0][i]
458  + stdVelocity[1][i]*stdVelocity[1][i]
459  + stdVelocity[2][i]*stdVelocity[2][i]);
460 
461  if (pntVelocity > stdV[el])
462  {
463  stdV[el] = pntVelocity;
464  }
465  }
466  }
467  }
468 
469  return stdV;
470  }
471 }
f SOLVER_UTILS_EXPORT void AdvectionNonConservativeForm(const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
Compute the non-conservative advection.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
static NekDouble EigenvaluesAnaMeshesRK2[6][14]
Definition: CFLtester.h:65
virtual void v_InitObject()
Init object for UnsteadySystem class.
Definition: CFLtester.cpp:54
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:848
int m_expdim
Expansion dimension.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Array< OneD, NekDouble > GetStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
Definition: CFLtester.cpp:329
virtual void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Definition: CFLtester.cpp:222
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: CellModel.h:51
virtual ~CFLtester()
Definition: CFLtester.cpp:78
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
SOLVER_UTILS_EXPORT int GetTotPoints()
static NekDouble EigenvaluesAnaMeshesRK4[6][14]
Definition: CFLtester.h:73
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
static std::string className
Definition: CFLtester.h:91
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition: CFLtester.cpp:82
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
boost::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:48
virtual NekDouble v_GetTimeStep(const Array< OneD, int > ExpOrder, const Array< OneD, NekDouble > CFL, NekDouble timeCFL)
Definition: CFLtester.cpp:229
int m_spacedim
Spatial dimension (>= expansion dim).
virtual void v_GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: CFLtester.cpp:172
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
EquationSystemFactory & GetEquationSystemFactory()
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition: CFLtester.cpp:131
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
Array< OneD, Array< OneD, NekDouble > > m_velocity
Definition: CFLtester.h:97
static NekDouble EigenvaluesAnaMeshesAB2[6][14]
Definition: CFLtester.h:57
CFLtester(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: CFLtester.cpp:48
SOLVER_UTILS_EXPORT void WeakDGAdvection(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
Calculate the weak discontinuous Galerkin advection.
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: CFLtester.h:86
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Geometry is curved or has non-constant factors.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
virtual void v_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
Definition: CFLtester.cpp:189
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215