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