Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UnsteadySystem.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File UnsteadySystem.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: Generic timestepping for Unsteady solvers
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <iomanip>
38 
43 
44 namespace Nektar
45 {
46  namespace SolverUtils
47  {
48  /**
49  * @class UnsteadySystem
50  *
51  * Provides the underlying timestepping framework for unsteady solvers
52  * including the general timestepping routines. This class is not
53  * intended to be directly instantiated, but rather is a base class
54  * on which to define unsteady solvers.
55  *
56  * For details on implementing unsteady solvers see
57  * \ref sectionADRSolverModuleImplementation here
58  */
59 
60  /**
61  * Processes SolverInfo parameters from the session file and sets up
62  * timestepping-specific code.
63  * @param pSession Session object to read parameters from.
64  */
67  : EquationSystem(pSession),
68  m_infosteps(10)
69 
70  {
71  }
72 
73  /**
74  * Initialization object for UnsteadySystem class.
75  */
77  {
79 
80  m_initialStep = 0;
81 
82  // Load SolverInfo parameters
83  m_session->MatchSolverInfo("DIFFUSIONADVANCEMENT","Explicit",
84  m_explicitDiffusion,true);
85  m_session->MatchSolverInfo("ADVECTIONADVANCEMENT","Explicit",
86  m_explicitAdvection,true);
87  m_session->MatchSolverInfo("REACTIONADVANCEMENT", "Explicit",
88  m_explicitReaction, true);
89 
90  // For steady problems, we do not initialise the time integration
91  if (m_session->DefinesSolverInfo("TIMEINTEGRATIONMETHOD"))
92  {
94  CreateInstance(m_session->GetSolverInfo(
95  "TIMEINTEGRATIONMETHOD"));
96 
97  // Load generic input parameters
98  m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
99  m_session->LoadParameter("CFL", m_cflSafetyFactor, 0.0);
100 
101  // Set up time to be dumped in field information
102  m_fieldMetaDataMap["Time"] =
103  boost::lexical_cast<std::string>(m_time);
104  }
105 
106  // By default attempt to forward transform initial condition.
107  m_homoInitialFwd = true;
108 
109  // Set up filters
110  LibUtilities::FilterMap::const_iterator x;
111  LibUtilities::FilterMap f = m_session->GetFilters();
112  for (x = f.begin(); x != f.end(); ++x)
113  {
114  m_filters.push_back(GetFilterFactory().CreateInstance(
115  x->first,
116  m_session,
117  x->second));
118  }
119  }
120 
121  /**
122  * Destructor for the class UnsteadyAdvection.
123  */
125  {
126  }
127 
128  /**
129  * @brief Returns the maximum time estimator for CFL control.
130  */
132  {
133  NekDouble TimeStability = 0.0;
134  switch(m_intScheme->GetIntegrationMethod())
135  {
139  {
140  TimeStability = 2.784;
141  break;
142  }
149  {
150  TimeStability = 2.0;
151  break;
152  }
154  {
155  TimeStability = 1.0;
156  break;
157  }
158  default:
159  {
160  ASSERTL0(
161  false,
162  "No CFL control implementation for this time"
163  "integration scheme");
164  }
165  }
166  return TimeStability;
167  }
168 
169  /**
170  * @brief Initialises the time integration scheme (as specified in the
171  * session file), and perform the time integration.
172  */
174  {
175  ASSERTL0(m_intScheme != 0, "No time integration scheme.");
176 
177  int i = 1;
178  int nvariables = 0;
179  int nfields = m_fields.num_elements();
180 
181  if (m_intVariables.empty())
182  {
183  for (i = 0; i < nfields; ++i)
184  {
185  m_intVariables.push_back(i);
186  }
187  nvariables = nfields;
188  }
189  else
190  {
191  nvariables = m_intVariables.size();
192  }
193 
194  // Integrate in wave-space if using homogeneous1D
196  {
197  for(i = 0; i < nfields; ++i)
198  {
199  m_fields[i]->HomogeneousFwdTrans(m_fields[i]->GetPhys(),
200  m_fields[i]->UpdatePhys());
201  m_fields[i]->SetWaveSpace(true);
202  m_fields[i]->SetPhysState(false);
203  }
204  }
205 
206  // Set up wrapper to fields data storage.
207  Array<OneD, Array<OneD, NekDouble> > fields(nvariables);
208  Array<OneD, Array<OneD, NekDouble> > tmp (nvariables);
209 
210  // Order storage to list time-integrated fields first.
211  for(i = 0; i < nvariables; ++i)
212  {
213  fields[i] = m_fields[m_intVariables[i]]->GetPhys();
214  m_fields[m_intVariables[i]]->SetPhysState(false);
215  }
216 
217  // Initialise time integration scheme
218  m_intSoln = m_intScheme->InitializeScheme(
219  m_timestep, fields, m_time, m_ode);
220 
221  // Initialise filters
223  for (x = m_filters.begin(); x != m_filters.end(); ++x)
224  {
225  (*x)->Initialise(m_fields, m_time);
226  }
227 
228  // Ensure that there is no conflict of parameters
229  if(m_cflSafetyFactor > 0.0)
230  {
231  // Check final condition
232  ASSERTL0(m_fintime == 0.0 || m_steps == 0,
233  "Final condition not unique: "
234  "fintime > 0.0 and Nsteps > 0");
235 
236  // Check timestep condition
237  ASSERTL0(m_timestep == 0.0,
238  "Timestep not unique: timestep > 0.0 & CFL > 0.0");
239  }
240 
241  // Check uniqueness of checkpoint output
242  ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
243  (m_checktime > 0.0 && m_checksteps == 0) ||
244  (m_checktime == 0.0 && m_checksteps > 0),
245  "Only one of IO_CheckTime and IO_CheckSteps "
246  "should be set!");
247 
248  Timer timer;
249  bool doCheckTime = false;
250  int step = m_initialStep;
251  NekDouble intTime = 0.0;
252  NekDouble lastCheckTime = 0.0;
253  NekDouble cpuTime = 0.0;
254  NekDouble elapsed = 0.0;
255 
256  while (step < m_steps ||
258  {
259  if (m_cflSafetyFactor)
260  {
261  m_timestep = GetTimeStep(fields);
262 
263  // Ensure that the final timestep finishes at the final
264  // time, or at a prescribed IO_CheckTime.
265  if (m_time + m_timestep > m_fintime && m_fintime > 0.0)
266  {
268  }
269  else if (m_checktime &&
270  m_time + m_timestep - lastCheckTime >= m_checktime)
271  {
272  lastCheckTime += m_checktime;
273  m_timestep = lastCheckTime - m_time;
274  doCheckTime = true;
275  }
276  }
277 
278  // Perform any solver-specific pre-integration steps
279  timer.Start();
280  if (v_PreIntegrate(step))
281  {
282  break;
283  }
284 
285  fields = m_intScheme->TimeIntegrate(
286  step, m_timestep, m_intSoln, m_ode);
287  timer.Stop();
288 
289  m_time += m_timestep;
290  elapsed = timer.TimePerTest(1);
291  intTime += elapsed;
292  cpuTime += elapsed;
293 
294  // Write out status information
295  if (m_session->GetComm()->GetRank() == 0 &&
296  !((step+1) % m_infosteps))
297  {
298  cout << "Steps: " << setw(8) << left << step+1 << " "
299  << "Time: " << setw(12) << left << m_time;
300 
301  if (m_cflSafetyFactor)
302  {
303  cout << " Time-step: " << setw(12)
304  << left << m_timestep;
305  }
306 
307  stringstream ss;
308  ss << cpuTime << "s";
309  cout << " CPU Time: " << setw(8) << left
310  << ss.str() << endl;
311  cpuTime = 0.0;
312  }
313 
314  // Transform data into coefficient space
315  for (i = 0; i < nvariables; ++i)
316  {
317  m_fields[m_intVariables[i]]->SetPhys(fields[i]);
318  m_fields[m_intVariables[i]]->FwdTrans_IterPerExp(
319  fields[i],
320  m_fields[m_intVariables[i]]->UpdateCoeffs());
321  m_fields[m_intVariables[i]]->SetPhysState(false);
322  }
323 
324  // Perform any solver-specific post-integration steps
325  if (v_PostIntegrate(step))
326  {
327  break;
328  }
329 
330  // search for NaN and quit if found
331  bool nanFound = false;
332  for (i = 0; i < nvariables; ++i)
333  {
334  if (Vmath::Nnan(fields[i].num_elements(), fields[i], 1) > 0)
335  {
336  cout << "NaN found in variable \""
337  << m_session->GetVariable(i)
338  << "\", terminating" << endl;
339  nanFound = true;
340  }
341  }
342 
343  if (nanFound)
344  {
345  break;
346  }
347 
348  // Update filters
350  for (x = m_filters.begin(); x != m_filters.end(); ++x)
351  {
352  (*x)->Update(m_fields, m_time);
353  }
354 
355  // Write out checkpoint files
356  if ((m_checksteps && step && !((step + 1) % m_checksteps)) ||
357  doCheckTime)
358  {
360  {
361  vector<bool> transformed(nfields, false);
362  for(i = 0; i < nfields; i++)
363  {
364  if (m_fields[i]->GetWaveSpace())
365  {
366  m_fields[i]->SetWaveSpace(false);
367  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
368  m_fields[i]->UpdatePhys());
369  m_fields[i]->SetPhysState(true);
370  transformed[i] = true;
371  }
372  }
374  for(i = 0; i < nfields; i++)
375  {
376  if (transformed[i])
377  {
378  m_fields[i]->SetWaveSpace(true);
379  m_fields[i]->HomogeneousFwdTrans(
380  m_fields[i]->GetPhys(),
381  m_fields[i]->UpdatePhys());
382  m_fields[i]->SetPhysState(false);
383  }
384  }
385  }
386  else
387  {
389  }
390  doCheckTime = false;
391  }
392 
393  // Step advance
394  ++step;
395  }
396 
397  // Print out summary statistics
398  if (m_session->GetComm()->GetRank() == 0)
399  {
400  if (m_cflSafetyFactor > 0.0)
401  {
402  cout << "CFL safety factor : " << m_cflSafetyFactor << endl
403  << "CFL time-step : " << m_timestep << endl;
404  }
405 
406  if (m_session->GetSolverInfo("Driver") != "SteadyState")
407  {
408  cout << "Time-integration : " << intTime << "s" << endl;
409  }
410  }
411 
412  // If homogeneous, transform back into physical space if necessary.
414  {
415  for(i = 0; i < nfields; i++)
416  {
417  if (m_fields[i]->GetWaveSpace())
418  {
419  m_fields[i]->SetWaveSpace(false);
420  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
421  m_fields[i]->UpdatePhys());
422  m_fields[i]->SetPhysState(true);
423  }
424  }
425  }
426  else
427  {
428  for(i = 0; i < nvariables; ++i)
429  {
430  m_fields[m_intVariables[i]]->SetPhys(fields[i]);
431  m_fields[m_intVariables[i]]->SetPhysState(true);
432  }
433  }
434 
435  // Finalise filters
436  for (x = m_filters.begin(); x != m_filters.end(); ++x)
437  {
438  (*x)->Finalise(m_fields, m_time);
439  }
440 
441  // Print for 1D problems
442  if(m_spacedim == 1)
443  {
444  v_AppendOutput1D(fields);
445  }
446  }
447 
448  /**
449  * @brief Sets the initial conditions.
450  */
452  {
456  }
457 
458  /**
459  * @brief Prints a summary with some information regards the
460  * time-stepping.
461  */
463  {
465  AddSummaryItem(s, "Advection",
466  m_explicitAdvection ? "explicit" : "implicit");
467 
468  if(m_session->DefinesSolverInfo("AdvectionType"))
469  {
470  AddSummaryItem(s, "AdvectionType",
471  m_session->GetSolverInfo("AdvectionType"));
472  }
473 
474  AddSummaryItem(s, "Diffusion",
475  m_explicitDiffusion ? "explicit" : "implicit");
476 
477  if (m_session->GetSolverInfo("EQTYPE")
478  == "SteadyAdvectionDiffusionReaction")
479  {
480  AddSummaryItem(s, "Reaction",
481  m_explicitReaction ? "explicit" : "implicit");
482  }
483 
484  AddSummaryItem(s, "Time Step", m_timestep);
485  AddSummaryItem(s, "No. of Steps", m_steps);
486  AddSummaryItem(s, "Checkpoints (steps)", m_checksteps);
487  AddSummaryItem(s, "Integration Type",
489  m_intScheme->GetIntegrationMethod()]);
490  }
491 
492  /**
493  * Stores the solution in a file for 1D problems only. This method has
494  * been implemented to facilitate the post-processing for 1D problems.
495  */
497  Array<OneD, Array<OneD, NekDouble> > &solution1D)
498  {
499  // Coordinates of the quadrature points in the real physical space
503  m_fields[0]->GetCoords(x, y, z);
504 
505  // Print out the solution in a txt file
506  ofstream outfile;
507  outfile.open("solution1D.txt");
508  for(int i = 0; i < GetNpoints(); i++)
509  {
510  outfile << scientific << setw (17) << setprecision(16) << x[i]
511  << " " << solution1D[0][i] << endl;
512  }
513  outfile << endl << endl;
514  outfile.close();
515  }
516 
518  Array<OneD, Array<OneD, NekDouble> > &physfield,
519  Array<OneD, Array<OneD, NekDouble> > &numflux)
520  {
521  ASSERTL0(false,
522  "This function is not implemented for this equation.");
523  }
524 
526  Array<OneD, Array<OneD, NekDouble> > &physfield,
527  Array<OneD, Array<OneD, NekDouble> > &numfluxX,
528  Array<OneD, Array<OneD, NekDouble> > &numfluxY )
529  {
530  ASSERTL0(false,
531  "This function is not implemented for this equation.");
532  }
533 
535  const Array<OneD, Array<OneD, NekDouble> > &ufield,
537  {
538  int i, j;
539  int nTraceNumPoints = GetTraceNpoints();
540  int nvariables = m_fields.num_elements();
541  int nqvar = uflux.num_elements();
542 
543  Array<OneD, NekDouble > Fwd (nTraceNumPoints);
544  Array<OneD, NekDouble > Bwd (nTraceNumPoints);
545  Array<OneD, NekDouble > Vn (nTraceNumPoints, 0.0);
546  Array<OneD, NekDouble > fluxtemp(nTraceNumPoints, 0.0);
547 
548  // Get the sign of (v \cdot n), v = an arbitrary vector
549 
550  // Evaulate upwind flux:
551  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
552  for (j = 0; j < nqvar; ++j)
553  {
554  for (i = 0; i < nvariables ; ++i)
555  {
556  // Compute Fwd and Bwd value of ufield of i direction
557  m_fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
558 
559  // if Vn >= 0, flux = uFwd, i.e.,
560  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
561  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
562 
563  // else if Vn < 0, flux = uBwd, i.e.,
564  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
565  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
566 
567  m_fields[i]->GetTrace()->Upwind(m_traceNormals[j],
568  Fwd, Bwd, fluxtemp);
569 
570  // Imposing weak boundary condition with flux
571  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
572  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
573  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
574 
575  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
576  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
577  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
578 
579  if(m_fields[0]->GetBndCondExpansions().num_elements())
580  {
581  WeakPenaltyforScalar(i, ufield[i], fluxtemp);
582  }
583 
584  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
585  // i.e,
586  // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
587  // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
588 
589  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
590  // i.e,
591  // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
592  // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
593 
594  Vmath::Vmul(nTraceNumPoints,
595  m_traceNormals[j], 1,
596  fluxtemp, 1,
597  uflux[j][i], 1);
598  }
599  }
600  }
601 
602 
603 
605  const Array<OneD, Array<OneD, NekDouble> > &ufield,
608  {
609  int nTraceNumPoints = GetTraceNpoints();
610  int nvariables = m_fields.num_elements();
611  int nqvar = qfield.num_elements();
612 
613  NekDouble C11 = 1.0;
614  Array<OneD, NekDouble > Fwd(nTraceNumPoints);
615  Array<OneD, NekDouble > Bwd(nTraceNumPoints);
616  Array<OneD, NekDouble > Vn (nTraceNumPoints, 0.0);
617 
618  Array<OneD, NekDouble > qFwd (nTraceNumPoints);
619  Array<OneD, NekDouble > qBwd (nTraceNumPoints);
620  Array<OneD, NekDouble > qfluxtemp(nTraceNumPoints, 0.0);
621 
622  Array<OneD, NekDouble > uterm(nTraceNumPoints);
623 
624  // Evaulate upwind flux:
625  // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
626  for (int i = 0; i < nvariables; ++i)
627  {
628  qflux[i] = Array<OneD, NekDouble> (nTraceNumPoints, 0.0);
629  for (int j = 0; j < nqvar; ++j)
630  {
631  // Compute Fwd and Bwd value of ufield of jth direction
632  m_fields[i]->GetFwdBwdTracePhys(qfield[j][i],qFwd,qBwd);
633 
634  // if Vn >= 0, flux = uFwd, i.e.,
635  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
636  // qflux = qBwd = q+
637  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
638  // qflux = qBwd = q-
639 
640  // else if Vn < 0, flux = uBwd, i.e.,
641  // edge::eForward, if V*n<0 <=> V*n_F<0, pick
642  // qflux = qFwd = q-
643  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
644  // qflux = qFwd = q+
645 
646  m_fields[i]->GetTrace()->Upwind(m_traceNormals[j],
647  qBwd, qFwd,
648  qfluxtemp);
649 
650  Vmath::Vmul(nTraceNumPoints,
651  m_traceNormals[j], 1,
652  qfluxtemp, 1,
653  qfluxtemp, 1);
654 
655  // Generate Stability term = - C11 ( u- - u+ )
656  m_fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
657 
658  Vmath::Vsub(nTraceNumPoints,
659  Fwd, 1, Bwd, 1,
660  uterm, 1);
661 
662  Vmath::Smul(nTraceNumPoints,
663  -1.0 * C11, uterm, 1,
664  uterm, 1);
665 
666  // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
667  Vmath::Vadd(nTraceNumPoints,
668  uterm, 1,
669  qfluxtemp, 1,
670  qfluxtemp, 1);
671 
672  // Imposing weak boundary condition with flux
673  if (m_fields[0]->GetBndCondExpansions().num_elements())
674  {
675  WeakPenaltyforVector(i, j,
676  qfield[j][i],
677  qfluxtemp,
678  C11);
679  }
680 
681  // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
682  // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
683  // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
684  Vmath::Vadd(nTraceNumPoints,
685  qfluxtemp, 1,
686  qflux[i], 1,
687  qflux[i], 1);
688  }
689  }
690  }
691 
693  {
694  if (m_session->DefinesFunction("InitialConditions"))
695  {
696  for (int i = 0; i < m_fields.num_elements(); ++i)
697  {
699 
700  vType = m_session->GetFunctionType(
701  "InitialConditions", m_session->GetVariable(i));
702 
703  if (vType == LibUtilities::eFunctionTypeFile)
704  {
705  std::string filename
706  = m_session->GetFunctionFilename(
707  "InitialConditions", m_session->GetVariable(i));
708 
709  fs::path pfilename(filename);
710 
711  // redefine path for parallel file which is in directory
712  if(fs::is_directory(pfilename))
713  {
714  fs::path metafile("Info.xml");
715  fs::path fullpath = pfilename / metafile;
716  filename = LibUtilities::PortablePath(fullpath);
717  }
718  m_fld->ImportFieldMetaData(filename, m_fieldMetaDataMap);
719 
720  // check to see if time defined
721  if (m_fieldMetaDataMap !=
723  {
725 
726  iter = m_fieldMetaDataMap.find("Time");
727  if (iter != m_fieldMetaDataMap.end())
728  {
729  time = boost::lexical_cast<NekDouble>(
730  iter->second);
731  }
732  }
733 
734  break;
735  }
736  }
737  }
738  }
739 
741  const int var,
742  const Array<OneD, const NekDouble> &physfield,
743  Array<OneD, NekDouble> &penaltyflux,
744  NekDouble time)
745  {
746  int i, e, npoints, id1, id2;
747 
748  // Number of boundary regions
749  int nbnd = m_fields[var]->GetBndCondExpansions().num_elements();
750  int Nfps, numBDEdge;
751  int nTraceNumPoints = GetTraceNpoints();
752  int cnt = 0;
753 
754  Array<OneD, NekDouble > uplus(nTraceNumPoints);
755 
756  m_fields[var]->ExtractTracePhys(physfield, uplus);
757  for (i = 0; i < nbnd; ++i)
758  {
759  // Number of boundary expansion related to that region
760  numBDEdge = m_fields[var]->
761  GetBndCondExpansions()[i]->GetExpSize();
762 
763  // Evaluate boundary values g_D or g_N from input files
765  m_session->GetFunction("InitialConditions", 0);
766 
767  npoints = m_fields[var]->
768  GetBndCondExpansions()[i]->GetNpoints();
769 
770  Array<OneD,NekDouble> BDphysics(npoints);
771  Array<OneD,NekDouble> x0(npoints,0.0);
772  Array<OneD,NekDouble> x1(npoints,0.0);
773  Array<OneD,NekDouble> x2(npoints,0.0);
774 
775  m_fields[var]->GetBndCondExpansions()[i]->GetCoords(x0,x1,x2);
776  ifunc->Evaluate(x0,x1,x2,time,BDphysics);
777 
778  // Weakly impose boundary conditions by modifying flux values
779  for (e = 0; e < numBDEdge ; ++e)
780  {
781  // Number of points on the expansion
782  Nfps = m_fields[var]->
783  GetBndCondExpansions()[i]->GetExp(e)->GetNumPoints(0);
784 
785  id1 = m_fields[var]->
786  GetBndCondExpansions()[i]->GetPhys_Offset(e);
787 
788  id2 = m_fields[0]->GetTrace()->
789  GetPhys_Offset(m_fields[0]->GetTraceMap()->
790  GetBndCondTraceToGlobalTraceMap(cnt++));
791 
792  // For Dirichlet boundary condition: uflux = g_D
793  if (m_fields[var]->GetBndConditions()[i]->
794  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
795  {
796  Vmath::Vcopy(Nfps,
797  &BDphysics[id1], 1,
798  &penaltyflux[id2], 1);
799  }
800  // For Neumann boundary condition: uflux = u+
801  else if ((m_fields[var]->GetBndConditions()[i])->
802  GetBoundaryConditionType() == SpatialDomains::eNeumann)
803  {
804  Vmath::Vcopy(Nfps,
805  &uplus[id2], 1,
806  &penaltyflux[id2], 1);
807  }
808  }
809  }
810  }
811 
812  /**
813  * Diffusion: Imposing weak boundary condition for q with flux
814  * uflux = g_D on Dirichlet boundary condition
815  * uflux = u_Fwd on Neumann boundary condition
816  */
818  const int var,
819  const int dir,
820  const Array<OneD, const NekDouble> &physfield,
821  Array<OneD, NekDouble> &penaltyflux,
822  NekDouble C11,
823  NekDouble time)
824  {
825  int i, e, npoints, id1, id2;
826  int nbnd = m_fields[var]->GetBndCondExpansions().num_elements();
827  int numBDEdge, Nfps;
828  int nTraceNumPoints = GetTraceNpoints();
829  Array<OneD, NekDouble > uterm(nTraceNumPoints);
830  Array<OneD, NekDouble > qtemp(nTraceNumPoints);
831  int cnt = 0;
832 
833  m_fields[var]->ExtractTracePhys(physfield,qtemp);
834 
835  for (i = 0; i < nbnd; ++i)
836  {
837  numBDEdge = m_fields[var]->
838  GetBndCondExpansions()[i]->GetExpSize();
839 
840  // Evaluate boundary values g_D or g_N from input files
842  m_session->GetFunction("InitialConditions", 0);
843 
844  npoints = m_fields[var]->
845  GetBndCondExpansions()[i]->GetNpoints();
846 
847  Array<OneD,NekDouble> BDphysics(npoints);
848  Array<OneD,NekDouble> x0(npoints,0.0);
849  Array<OneD,NekDouble> x1(npoints,0.0);
850  Array<OneD,NekDouble> x2(npoints,0.0);
851 
852  m_fields[var]->GetBndCondExpansions()[i]->GetCoords(x0,x1,x2);
853  ifunc->Evaluate(x0,x1,x2,time,BDphysics);
854 
855  // Weakly impose boundary conditions by modifying flux values
856  for (e = 0; e < numBDEdge ; ++e)
857  {
858  Nfps = m_fields[var]->
859  GetBndCondExpansions()[i]->GetExp(e)->GetNumPoints(0);
860 
861  id1 = m_fields[var]->
862  GetBndCondExpansions()[i]->GetPhys_Offset(e);
863 
864  id2 = m_fields[0]->GetTrace()->
865  GetPhys_Offset(m_fields[0]->GetTraceMap()->
866  GetBndCondTraceToGlobalTraceMap(cnt++));
867 
868  // For Dirichlet boundary condition:
869  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
870  if(m_fields[var]->GetBndConditions()[i]->
871  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
872  {
873  Vmath::Vmul(Nfps,
874  &m_traceNormals[dir][id2], 1,
875  &qtemp[id2], 1,
876  &penaltyflux[id2], 1);
877  }
878  // For Neumann boundary condition: qflux = g_N
879  else if((m_fields[var]->GetBndConditions()[i])->
880  GetBoundaryConditionType() == SpatialDomains::eNeumann)
881  {
882  Vmath::Vmul(Nfps,
883  &m_traceNormals[dir][id2], 1,
884  &BDphysics[id1], 1,
885  &penaltyflux[id2], 1);
886  }
887  }
888  }
889  }
890 
891  /**
892  * @brief Return the timestep to be used for the next step in the
893  * time-marching loop.
894  *
895  * This function can be overloaded to facilitate solver which utilise a
896  * CFL (or other) parameter to determine a maximum timestep under which
897  * the problem remains stable.
898  */
900  const Array<OneD, const Array<OneD, NekDouble> > &inarray)
901  {
902  return v_GetTimeStep(inarray);
903  }
904 
905  /**
906  * @brief Return the timestep to be used for the next step in the
907  * time-marching loop.
908  *
909  * @see UnsteadySystem::GetTimeStep
910  */
912  const Array<OneD, const Array<OneD, NekDouble> > &inarray)
913  {
914  ASSERTL0(false, "Not defined for this class");
915  return 0.0;
916  }
917 
919  {
920  return false;
921  }
922 
924  {
925  return false;
926  }
927 
929  {
930  return false;
931  }
932 
934  const Array<OneD, Array<OneD, NekDouble> > vel,
935  StdRegions::VarCoeffMap &varCoeffMap)
936  {
937  int phystot = m_fields[0]->GetTotPoints();
938  int nvel = vel.num_elements();
939 
940  Array<OneD, NekDouble> varcoeff(phystot),tmp;
941 
942  // calculate magnitude of v
943  Vmath::Vmul(phystot,vel[0],1,vel[0],1,varcoeff,1);
944  for(int n = 1; n < nvel; ++n)
945  {
946  Vmath::Vvtvp(phystot,vel[n],1,vel[n],1,varcoeff,1,varcoeff,1);
947  }
948  Vmath::Vsqrt(phystot,varcoeff,1,varcoeff,1);
949 
950  for(int i = 0; i < m_fields[0]->GetNumElmts(); ++i)
951  {
952  int offset = m_fields[0]->GetPhys_Offset(i);
953  int nq = m_fields[0]->GetExp(i)->GetTotPoints();
954  Array<OneD, NekDouble> unit(nq,1.0);
955 
956  int nmodes = 0;
957 
958  for(int n = 0; n < m_fields[0]->GetExp(i)->GetNumBases(); ++n)
959  {
960  nmodes = max(nmodes,
961  m_fields[0]->GetExp(i)->GetBasisNumModes(n));
962  }
963 
964  NekDouble h = m_fields[0]->GetExp(i)->Integral(unit);
965  h = pow(h,(NekDouble) (1.0/nvel))/((NekDouble) nmodes);
966 
967  Vmath::Smul(nq,h,varcoeff+offset,1,tmp = varcoeff+offset,1);
968  }
969 
970  // set up map with eVarCoffLaplacian key
971  varCoeffMap[StdRegions::eVarCoeffLaplacian] = varcoeff;
972  }
973  }
974 }
virtual SOLVER_UTILS_EXPORT bool v_SteadyStateCheck(int step)
virtual SOLVER_UTILS_EXPORT void v_DoSolve()
Solves an unsteady problem.
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &l)
Virtual function for generating summary information.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
bool m_explicitReaction
Indicates if explicit or implicit treatment of reaction is used.
virtual SOLVER_UTILS_EXPORT void v_NumFluxforVector(const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
SOLVER_UTILS_EXPORT void CheckForRestartTime(NekDouble &time)
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate(int step)
A base class for describing how to solve specific equations.
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
Adams-Bashforth Forward multi-step scheme of order 2.
NekDouble m_time
Current time of simulation.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Runge-Kutta multi-stage scheme 4th order explicit (old name)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_timestep
Time step size.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
void WeakPenaltyforVector(const int var, const int dir, const Array< OneD, const NekDouble > &physfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11, NekDouble time=0.0)
Classical RungeKutta4 method (new name for eClassicalRungeKutta4)
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
virtual SOLVER_UTILS_EXPORT void v_AppendOutput1D(Array< OneD, Array< OneD, NekDouble > > &solution1D)
Print the solution at each solution point in a txt file.
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
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep(const Array< OneD, const Array< OneD, NekDouble > > &inarray)
Return the timestep to be used for the next step in the time-marching loop.
NekDouble m_checktime
Time between checkpoints.
const char *const TimeIntegrationMethodMap[]
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator()
Get the maximum timestep estimator for cfl control.
Nonlinear SSP RungeKutta3 explicit.
int m_nchk
Number of checkpoints written so far.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
int m_checksteps
Number of steps between checkpoints.
NekDouble m_fintime
Finish time of the simulation.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Initialisation object for EquationSystem.
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff(const Array< OneD, Array< OneD, NekDouble > > vel, StdRegions::VarCoeffMap &varCoeffMap)
Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h t...
int m_steps
Number of steps to take.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Classical RungeKutta2 method (new name for eMidpoint)
void Stop()
Definition: Timer.cpp:62
static const NekDouble kNekZeroTol
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.
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem()
Destructor.
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate(int step)
Adams-Bashforth Forward multi-step scheme of order 1.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:226
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 Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.cpp:878
std::vector< std::pair< std::string, FilterParams > > FilterMap
Definition: SessionReader.h:66
Nonlinear SSP RungeKutta2 explicit (surrogate for eRungeKutta2_ImprovedEuler)
Improved RungeKutta2 explicit (old name meaning Heun's method)
int m_spacedim
Spatial dimension (>= expansion dim).
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
std::string PortablePath(const boost::filesystem::path &path)
create portable path on different platforms for boost::filesystem path
Definition: FileSystem.cpp:41
double NekDouble
TimeIntegrationWrapperFactory & GetTimeIntegrationWrapperFactory()
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT NekDouble GetTimeStep()
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
virtual SOLVER_UTILS_EXPORT void v_NumFluxforScalar(const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
boost::shared_ptr< Equation > EquationSharedPtr
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession)
Initialises UnsteadySystem class members.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
virtual SOLVER_UTILS_EXPORT void v_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
Wrapper to the time integration scheme.
LibUtilities::FieldIOSharedPtr m_fld
Field input/output.
int m_initialStep
Number of the step where the simulation should begin.
FilterFactory & GetFilterFactory()
Definition: Filter.cpp:42
std::vector< FilterSharedPtr > m_filters
int m_infosteps
Number of time steps between outputting status information.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
void Start()
Definition: Timer.cpp:51
void WeakPenaltyforScalar(const int var, const Array< OneD, const NekDouble > &physfield, Array< OneD, NekDouble > &penaltyflux, NekDouble time=0.0)
midpoint method (old name)
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:108
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:54
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Sets up initial conditions.
enum HomogeneousType m_HomogeneousType
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
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