Nektar++
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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Generic timestepping for Unsteady solvers
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36#include <iostream>
37using namespace std;
38
39#include <boost/format.hpp>
40
45
46namespace Nektar::SolverUtils
47{
50 "set-start-time", "", "Set the starting time of the simulation.");
51
54 "set-start-chknumber", "",
55 "Set the starting number of the checkpoint file.");
56
57/**
58 * @class UnsteadySystem
59 *
60 * Provides the underlying timestepping framework for unsteady solvers
61 * including the general timestepping routines. This class is not
62 * intended to be directly instantiated, but rather is a base class
63 * on which to define unsteady solvers.
64 *
65 * For details on implementing unsteady solvers see
66 * \ref sectionADRSolverModuleImplementation here
67 */
68
69/**
70 * Processes SolverInfo parameters from the session file and sets up
71 * timestepping-specific code.
72 * @param pSession Session object to read parameters from.
73 */
77 : EquationSystem(pSession, pGraph), SolverUtils::ALEHelper()
78
79{
80}
81
82/**
83 * Initialization object for UnsteadySystem class.
84 */
85void UnsteadySystem::v_InitObject(bool DeclareField)
86{
87 EquationSystem::v_InitObject(DeclareField);
89 m_initialStep = 0;
90
91 // Load SolverInfo parameters.
92 m_session->MatchSolverInfo("DIFFUSIONADVANCEMENT", "Explicit",
94 m_session->MatchSolverInfo("ADVECTIONADVANCEMENT", "Explicit",
96 m_session->MatchSolverInfo("REACTIONADVANCEMENT", "Explicit",
97 m_explicitReaction, true);
98 m_session->LoadParameter("CheckAbortSteps", m_abortSteps, 1);
99
100 // Steady state tolerance.
101 m_session->LoadParameter("SteadyStateTol", m_steadyStateTol, 0.0);
102
103 // Frequency for checking steady state.
104 m_session->LoadParameter("SteadyStateSteps", m_steadyStateSteps, 1);
105
106 // For steady problems, we do not initialise the time integration.
107 if (m_session->DefinesSolverInfo("TimeIntegrationMethod") ||
108 m_session->DefinesTimeIntScheme())
109 {
111 if (m_session->DefinesTimeIntScheme())
112 {
113 timeInt = m_session->GetTimeIntScheme();
114 }
115 else
116 {
117 std::cout << "TimeIntegrationMethod is deprecated, please use "
118 "TIMEINTEGRATIONSCHEME";
119 timeInt.method = m_session->GetSolverInfo("TimeIntegrationMethod");
120 }
121
124 timeInt.method, timeInt.variant, timeInt.order,
125 timeInt.freeParams);
126
127 // Load generic input parameters.
128 m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
129 m_session->LoadParameter("IO_FiltersInfoSteps", m_filtersInfosteps,
130 10 * m_infosteps);
131 m_session->LoadParameter("CFL", m_cflSafetyFactor, 0.0);
132 m_session->LoadParameter("CFLEnd", m_CFLEnd, 0.0);
133 m_session->LoadParameter("CFLGrowth", m_CFLGrowth, 1.0);
134
135 // Ensure that there is no conflict of parameters.
136 if (m_cflSafetyFactor > 0.0)
137 {
138 // Check final condition.
139 ASSERTL0(m_fintime == 0.0 || m_steps == 0,
140 "Final condition not unique: "
141 "fintime > 0.0 and Nsteps > 0");
142 // Check timestep condition.
143 ASSERTL0(m_timestep == 0.0,
144 "Timestep not unique: timestep > 0.0 & CFL > 0.0");
145 }
146 else
147 {
148 ASSERTL0(m_timestep != 0.0, "Need to set either TimeStep or CFL");
149 }
150
151 // Ensure that there is no conflict of parameters.
152 if (m_CFLGrowth > 1.0)
153 {
154 // Check final condition.
156 "m_CFLEnd >= m_cflSafetyFactor required");
157 }
158
159 // Set up time to be dumped in field information.
160 m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
161 }
162
163 // By default attempt to forward transform initial condition.
164 m_homoInitialFwd = true;
165
166 // Set up filters.
167 for (auto &x : m_session->GetFilters())
168 {
169 m_filters.push_back(make_pair(
170 x.name, GetFilterFactory().CreateInstance(
171 x.name, m_session, shared_from_this(), x.params)));
172 }
173}
174
175/**
176 * @brief Returns the maximum time estimator for CFL control.
177 */
179{
180 return m_intScheme->GetTimeStability();
181}
182
183/**
184 * @brief Returns the time integration scheme.
185 */
188{
189 return m_intScheme;
190}
191
192/**
193 * @brief Returns the time integration scheme operators.
194 */
197{
198 return m_ode;
199}
200
201/**
202 * @brief Initialises the time integration scheme (as specified in the
203 * session file), and perform the time integration.
204 */
206{
207 ASSERTL0(m_intScheme != nullptr, "No time integration scheme.");
208
209 int i = 1;
210 int nvariables = 0;
211 int nfields = m_fields.size();
212 if (m_intVariables.empty())
213 {
214 for (i = 0; i < nfields; ++i)
215 {
216 m_intVariables.push_back(i);
217 }
218 nvariables = nfields;
219 }
220 else
221 {
222 nvariables = m_intVariables.size();
223 }
224
225 // Integrate in wave-space if using homogeneous1D.
227 {
228 for (i = 0; i < nfields; ++i)
229 {
230 m_fields[i]->HomogeneousFwdTrans(m_fields[i]->GetTotPoints(),
231 m_fields[i]->GetPhys(),
232 m_fields[i]->UpdatePhys());
233 m_fields[i]->SetWaveSpace(true);
234 m_fields[i]->SetPhysState(false);
235 }
236 }
237
238 // Set up wrapper to fields data storage.
239 Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
240
241 // Order storage to list time-integrated fields first.
242 // @TODO: Refactor to take coeffs (FwdTrans) if boolean flag (in constructor
243 // function) says to.
244 for (i = 0; i < nvariables; ++i)
245 {
246 fields[i] = m_fields[m_intVariables[i]]->UpdatePhys();
247 m_fields[m_intVariables[i]]->SetPhysState(false);
248 }
249
250 // @TODO: Virtual function that allows to transform the field space, embed
251 // the MultiplyMassMatrix in here.
252 // @TODO: Specify what the fields variables are physical or coefficient,
253 // boolean in UnsteadySystem class...
254
255 v_ALEPreMultiplyMass(fields);
256
257 // Initialise time integration scheme.
258 m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
259
260 // Initialise filters.
261 for (auto &x : m_filters)
262 {
263 x.second->Initialise(m_fields, m_time);
264 }
265
267 bool doCheckTime = false;
268 int step = m_initialStep;
269 int stepCounter = 0;
270 NekDouble intTime = 0.0;
271 NekDouble cpuTime = 0.0;
272 NekDouble cpuPrevious = 0.0;
273 NekDouble elapsed = 0.0;
274 NekDouble totFilterTime = 0.0;
275
276 m_lastCheckTime = 0.0;
277
278 Array<OneD, int> abortFlags(2, 0);
279 string abortFile = "abort";
280 if (m_session->DefinesSolverInfo("CheckAbortFile"))
281 {
282 abortFile = m_session->GetSolverInfo("CheckAbortFile");
283 }
284
285 NekDouble tmp_cflSafetyFactor = m_cflSafetyFactor;
286 while ((step < m_steps || m_time < m_fintime - NekConstants::kNekZeroTol) &&
287 abortFlags[1] == 0)
288 {
290 {
291 tmp_cflSafetyFactor =
292 min(m_CFLEnd, m_CFLGrowth * tmp_cflSafetyFactor);
293 }
294
295 // Frozen preconditioner checks.
296 if (!m_comm->IsParallelInTime())
297 {
299 {
300 m_cflSafetyFactor = tmp_cflSafetyFactor;
301
303 {
304 m_timestep = GetTimeStep(fields);
305 }
306
307 // Ensure that the final timestep finishes at the final
308 // time, or at a prescribed IO_CheckTime.
309 if (m_time + m_timestep > m_fintime && m_fintime > 0.0)
310 {
312 }
313 else if (m_checktime &&
315 {
318 doCheckTime = true;
319 }
320 }
321 }
322
323 if (m_TimeIncrementFactor > 1.0)
324 {
325 NekDouble timeincrementFactor = m_TimeIncrementFactor;
326 m_timestep *= timeincrementFactor;
327
328 if (m_time + m_timestep > m_fintime && m_fintime > 0.0)
329 {
331 }
332 }
333
334 // Perform any solver-specific pre-integration steps.
335 timer.Start();
336 if (v_PreIntegrate(
337 step)) // Could be possible to put a preintegrate step in the
338 // ALEHelper class, put in the Unsteady Advection class
339 {
340 break;
341 }
342
343 ASSERTL0(m_timestep > 0, "m_timestep < 0");
344
345 fields = m_intScheme->TimeIntegrate(stepCounter, m_timestep);
346 timer.Stop();
347
349 elapsed = timer.TimePerTest(1);
350 intTime += elapsed;
351 cpuTime += elapsed;
352
353 // Write out status information.
354 v_PrintStatusInformation(step, cpuTime);
355 if (m_infosteps &&
356 m_session->GetComm()->GetSpaceComm()->GetRank() == 0 &&
357 !((step + 1) % m_infosteps))
358 {
359 cpuPrevious = cpuTime;
360 cpuTime = 0.0;
361 }
362
363 // @TODO: Another virtual function with this in it based on if ALE or
364 // not.
365 if (m_ALESolver) // Change to advect coeffs, change flag to physical vs
366 // coefficent space
367 {
370 }
371 else
372 {
373 // Transform data into coefficient space
374 for (i = 0; i < nvariables; ++i)
375 {
376 // copy fields into ExpList::m_phys and assign the new
377 // array to fields
378 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
379 fields[i] = m_fields[m_intVariables[i]]->UpdatePhys();
380 if (v_RequireFwdTrans())
381 {
382 if (m_comm->IsParallelInTime())
383 {
384 m_fields[m_intVariables[i]]->FwdTrans(
385 m_fields[m_intVariables[i]]->GetPhys(),
386 m_fields[m_intVariables[i]]->UpdateCoeffs());
387 }
388 else
389 {
390 m_fields[m_intVariables[i]]->FwdTransLocalElmt(
391 m_fields[m_intVariables[i]]->GetPhys(),
392 m_fields[m_intVariables[i]]->UpdateCoeffs());
393 }
394 }
395 m_fields[m_intVariables[i]]->SetPhysState(false);
396 }
397 }
398
399 // Perform any solver-specific post-integration steps.
400 if (v_PostIntegrate(step))
401 {
402 break;
403 }
404
405 // Check for steady-state.
407 (!((step + 1) % m_steadyStateSteps)))
408 {
409 if (CheckSteadyState(step, intTime))
410 {
411 if (m_comm->GetRank() == 0)
412 {
413 cout << "Reached Steady State to tolerance "
414 << m_steadyStateTol << endl;
415 }
416 break;
417 }
418 }
419
420 // Test for abort conditions (nan, or abort file).
421 if (m_abortSteps && !((step + 1) % m_abortSteps))
422 {
423 abortFlags[0] = 0;
424 for (i = 0; i < nvariables; ++i)
425 {
426 if (Vmath::Nnan(m_fields[m_intVariables[i]]->GetPhys().size(),
427 m_fields[m_intVariables[i]]->GetPhys(), 1) > 0)
428 {
429 abortFlags[0] = 1;
430 }
431 }
432
433 // Rank zero looks for abort file and deletes it
434 // if it exists. Communicates the abort.
435 if (m_session->GetComm()->GetSpaceComm()->GetRank() == 0)
436 {
437 if (fs::exists(abortFile))
438 {
439 fs::remove(abortFile);
440 abortFlags[1] = 1;
441 }
442 }
443
444 m_session->GetComm()->GetSpaceComm()->AllReduce(
445 abortFlags, LibUtilities::ReduceMax);
446
447 ASSERTL0(!abortFlags[0], "NaN found during time integration.");
448 }
449
450 // Update filters.
451 for (auto &x : m_filters)
452 {
453 timer.Start();
454 x.second->Update(m_fields, m_time);
455 timer.Stop();
456 elapsed = timer.TimePerTest(1);
457 totFilterTime += elapsed;
458
459 // Write out individual filter status information.
460 if (m_filtersInfosteps && m_session->GetComm()->GetRank() == 0 &&
461 !((step + 1) % m_filtersInfosteps) && !m_filters.empty() &&
462 m_session->DefinesCmdLineArgument("verbose"))
463 {
464 stringstream s0;
465 s0 << x.first << ":";
466 stringstream s1;
467 s1 << elapsed << "s";
468 stringstream s2;
469 s2 << elapsed / cpuPrevious * 100 << "%";
470 cout << "CPU time for filter " << setw(25) << left << s0.str()
471 << setw(12) << left << s1.str() << endl
472 << "\t Percentage of time integration: " << setw(10)
473 << left << s2.str() << endl;
474 }
475 }
476
477 // Write out overall filter status information.
478 if (m_filtersInfosteps && m_session->GetComm()->GetRank() == 0 &&
479 !((step + 1) % m_filtersInfosteps) && !m_filters.empty())
480 {
481 stringstream ss;
482 ss << totFilterTime << "s";
483 cout << "Total filters CPU Time:\t\t\t " << setw(10) << left
484 << ss.str() << endl;
485 }
486 totFilterTime = 0.0;
487
488 // Write out checkpoint files.
489 if ((m_checksteps && !((step + 1) % m_checksteps)) || doCheckTime)
490 {
492 {
493 // Transform to physical space for output.
494 vector<bool> transformed(nfields, false);
495 for (i = 0; i < nfields; i++)
496 {
497 if (m_fields[i]->GetWaveSpace())
498 {
499 m_fields[i]->SetWaveSpace(false);
500 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
501 m_fields[i]->UpdatePhys());
502 transformed[i] = true;
503 }
504 }
506 m_nchk++;
507
508 // Transform back to wave-space after output.
509 for (i = 0; i < nfields; i++)
510 {
511 if (transformed[i])
512 {
513 m_fields[i]->SetWaveSpace(true);
514 m_fields[i]->HomogeneousFwdTrans(
515 m_fields[i]->GetTotPoints(), m_fields[i]->GetPhys(),
516 m_fields[i]->UpdatePhys());
517 m_fields[i]->SetPhysState(false);
518 }
519 }
520 }
521 else
522 {
524 m_nchk++;
525 }
526 doCheckTime = false;
527 }
528
529 // Step advance.
530 ++step;
531 ++stepCounter;
532 }
533
534 // Print out summary statistics.
536
537 // If homogeneous, transform back into physical space if necessary.
538 if (!m_ALESolver)
539 {
541 {
542 for (i = 0; i < nfields; i++)
543 {
544 if (m_fields[i]->GetWaveSpace())
545 {
546 m_fields[i]->SetWaveSpace(false);
547 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
548 m_fields[i]->UpdatePhys());
549 }
550 }
551 }
552 else
553 {
554 for (i = 0; i < nvariables; ++i)
555 {
556 m_fields[m_intVariables[i]]->SetPhysState(true);
557 }
558 }
559 }
560 // Finalise filters.
561 for (auto &x : m_filters)
562 {
563 x.second->Finalise(m_fields, m_time);
564 }
565
566 // Print for 1D problems.
567 if (m_spacedim == 1)
568 {
570 }
571}
572
574 const NekDouble cpuTime)
575{
576 if (m_infosteps && m_session->GetComm()->GetSpaceComm()->GetRank() == 0 &&
577 !((step + 1) % m_infosteps))
578 {
579 if (m_comm->IsParallelInTime())
580 {
581 cout << "RANK " << m_session->GetComm()->GetTimeComm()->GetRank()
582 << " Steps: " << setw(8) << left << step + 1 << " "
583 << "Time: " << setw(12) << left << m_time;
584 }
585 else
586 {
587 cout << "Steps: " << setw(8) << left << step + 1 << " "
588 << "Time: " << setw(12) << left << m_time;
589 }
590
592 {
593 cout << " Time-step: " << setw(12) << left << m_timestep;
594 }
595
596 stringstream ss;
597 ss << cpuTime << "s";
598 cout << " CPU Time: " << setw(8) << left << ss.str() << endl;
599 }
600}
601
603{
604 if (m_session->GetComm()->GetRank() == 0)
605 {
606 if (m_cflSafetyFactor > 0.0)
607 {
608 cout << "CFL safety factor : " << m_cflSafetyFactor << endl
609 << "CFL time-step : " << m_timestep << endl;
610 }
611
612 if (m_session->GetSolverInfo("Driver") != "SteadyState" &&
613 m_session->GetSolverInfo("Driver") != "Parareal" &&
614 m_session->GetSolverInfo("Driver") != "PFASST")
615 {
616 cout << "Time-integration : " << intTime << "s" << endl;
617 }
618 }
619}
620
621/**
622 * @brief Sets the initial conditions.
623 */
624void UnsteadySystem::v_DoInitialise(bool dumpInitialConditions)
625{
628 SetInitialConditions(m_time, dumpInitialConditions);
629
631
633}
634
635/**
636 * @brief Prints a summary with some information regards the
637 * time-stepping.
638 */
640{
642 AddSummaryItem(s, "Advect. advancement",
643 m_explicitAdvection ? "explicit" : "implicit");
644
645 AddSummaryItem(s, "Diffuse. advancement",
646 m_explicitDiffusion ? "explicit" : "implicit");
647
648 if (m_session->GetSolverInfo("EQTYPE") ==
649 "SteadyAdvectionDiffusionReaction")
650 {
651 AddSummaryItem(s, "React. advancement",
652 m_explicitReaction ? "explicit" : "implicit");
653 }
654
655 AddSummaryItem(s, "Time Step", m_timestep);
656 AddSummaryItem(s, "No. of Steps", m_steps);
657 AddSummaryItem(s, "Checkpoints (steps)", m_checksteps);
658 if (m_intScheme)
659 {
660 AddSummaryItem(s, "Integration Type", m_intScheme->GetName());
661 }
662}
663
664/**
665 * Stores the solution in a file for 1D problems only. This method has
666 * been implemented to facilitate the post-processing for 1D problems.
667 */
669{
670 // Coordinates of the quadrature points in the real physical space.
674 m_fields[0]->GetCoords(x, y, z);
675
676 // Print out the solution in a txt file.
677 ofstream outfile;
678 outfile.open("solution1D.txt");
679 for (int i = 0; i < GetNpoints(); i++)
680 {
681 outfile << scientific << setw(17) << setprecision(16) << x[i] << " "
682 << m_fields[m_intVariables[0]]->GetPhys()[i] << endl;
683 }
684 outfile << endl << endl;
685 outfile.close();
686}
687
688/**
689 *
690 */
692{
693 if (m_session->DefinesFunction("InitialConditions"))
694 {
695 for (int i = 0; i < m_fields.size(); ++i)
696 {
698
699 vType = m_session->GetFunctionType("InitialConditions",
700 m_session->GetVariable(i));
701
703 {
704 std::string filename = 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 }
718 fld->ImportFieldMetaData(filename, m_fieldMetaDataMap);
719
720 // Check to see if time defined.
722 {
723 auto iter = m_fieldMetaDataMap.find("Time");
724 if (iter != m_fieldMetaDataMap.end())
725 {
726 time = std::stod(iter->second);
727 }
728
729 iter = m_fieldMetaDataMap.find("ChkFileNum");
730 if (iter != m_fieldMetaDataMap.end())
731 {
732 nchk = std::stod(iter->second);
733 }
734 }
735
736 break;
737 }
738 }
739 }
740 if (m_session->DefinesCmdLineArgument("set-start-time"))
741 {
742 time = std::stod(
743 m_session->GetCmdLineArgument<std::string>("set-start-time")
744 .c_str());
745 }
746 if (m_session->DefinesCmdLineArgument("set-start-chknumber"))
747 {
748 nchk = std::stoi(
749 m_session->GetCmdLineArgument<std::string>("set-start-chknumber"));
750 }
751 ASSERTL0(time >= 0 && nchk >= 0,
752 "Starting time and checkpoint number should be >= 0");
753}
754
755/**
756 * @brief Return the timestep to be used for the next step in the
757 * time-marching loop.
758 *
759 * @see UnsteadySystem::GetTimeStep
760 */
762 [[maybe_unused]] const Array<OneD, const Array<OneD, NekDouble>> &inarray)
763{
764 NEKERROR(ErrorUtil::efatal, "Not defined for this class");
765 return 0.0;
766}
767
768/**
769 *
770 */
771bool UnsteadySystem::v_PreIntegrate([[maybe_unused]] int step)
772{
773 return false;
774}
775
776/**
777 *
778 */
779bool UnsteadySystem::v_PostIntegrate([[maybe_unused]] int step)
780{
781 return false;
782}
783
784/**
785 *
786 */
788{
789 return true;
790}
791
792/**
793 *
794 */
796{
797 return true;
798}
799
800/**
801 *
802 */
805 StdRegions::VarCoeffMap &varCoeffMap)
806{
807 int phystot = m_fields[0]->GetTotPoints();
808 int nvel = vel.size();
809
810 Array<OneD, NekDouble> varcoeff(phystot), tmp;
811
812 // Calculate magnitude of v.
813 Vmath::Vmul(phystot, vel[0], 1, vel[0], 1, varcoeff, 1);
814 for (int n = 1; n < nvel; ++n)
815 {
816 Vmath::Vvtvp(phystot, vel[n], 1, vel[n], 1, varcoeff, 1, varcoeff, 1);
817 }
818 Vmath::Vsqrt(phystot, varcoeff, 1, varcoeff, 1);
819
820 for (int i = 0; i < m_fields[0]->GetNumElmts(); ++i)
821 {
822 int offset = m_fields[0]->GetPhys_Offset(i);
823 int nq = m_fields[0]->GetExp(i)->GetTotPoints();
824 Array<OneD, NekDouble> unit(nq, 1.0);
825
826 int nmodes = 0;
827
828 for (int n = 0; n < m_fields[0]->GetExp(i)->GetNumBases(); ++n)
829 {
830 nmodes = max(nmodes, m_fields[0]->GetExp(i)->GetBasisNumModes(n));
831 }
832
833 NekDouble h = m_fields[0]->GetExp(i)->Integral(unit);
834 h = pow(h, (NekDouble)(1.0 / nvel)) / ((NekDouble)nmodes);
835
836 Vmath::Smul(nq, h, varcoeff + offset, 1, tmp = varcoeff + offset, 1);
837 }
838
839 // Set up map with eVarCoffLaplacian key.
840 varCoeffMap[StdRegions::eVarCoeffLaplacian] = varcoeff;
841}
842
843/**
844 *
845 */
847{
848 if (m_steadyStateTol > 0.0)
849 {
850 const int nPoints = m_fields[0]->GetTotPoints();
853
854 for (int i = 0; i < m_fields.size(); ++i)
855 {
857 Vmath::Vcopy(nPoints, m_fields[i]->GetPhys(), 1,
858 m_previousSolution[i], 1);
859 }
860
861 if (m_comm->GetRank() == 0)
862 {
863 std::string fName =
864 m_session->GetSessionName() + std::string(".resd");
865 m_errFile.open(fName.c_str());
866 m_errFile << setw(26) << left << "# Time";
867
868 m_errFile << setw(26) << left << "CPU_Time";
869
870 m_errFile << setw(26) << left << "Step";
871
872 for (int i = 0; i < m_fields.size(); ++i)
873 {
874 m_errFile << setw(26) << m_session->GetVariables()[i];
875 }
876
877 m_errFile << endl;
878 }
879 }
880}
881
882/**
883 * @brief Calculate whether the system has reached a steady state by
884 * observing residuals to a user-defined tolerance.
885 */
886bool UnsteadySystem::CheckSteadyState(int step, const NekDouble &totCPUTime)
887{
888 const int nPoints = GetTotPoints();
889 const int nFields = m_fields.size();
890
891 // Holds L2 errors.
892 Array<OneD, NekDouble> L2(nFields);
893
894 SteadyStateResidual(step, L2);
895
896 if (m_infosteps && m_comm->GetRank() == 0 &&
897 (((step + 1) % m_infosteps == 0) || ((step == m_initialStep))))
898 {
899 // Output time.
900 m_errFile << boost::format("%25.19e") % m_time;
901
902 m_errFile << " " << boost::format("%25.19e") % totCPUTime;
903
904 int stepp = step + 1;
905
906 m_errFile << " " << boost::format("%25.19e") % stepp;
907
908 // Output residuals.
909 for (int i = 0; i < nFields; ++i)
910 {
911 m_errFile << " " << boost::format("%25.19e") % L2[i];
912 }
913
914 m_errFile << endl;
915 }
916
917 // Calculate maximum L2 error.
918 NekDouble maxL2 = Vmath::Vmax(nFields, L2, 1);
919
920 if (m_infosteps && m_session->DefinesCmdLineArgument("verbose") &&
921 m_comm->GetRank() == 0 && ((step + 1) % m_infosteps == 0))
922 {
923 cout << "-- Maximum L^2 residual: " << maxL2 << endl;
924 }
925
926 if (maxL2 <= m_steadyStateTol)
927 {
928 return true;
929 }
930
931 for (int i = 0; i < m_fields.size(); ++i)
932 {
933 Vmath::Vcopy(nPoints, m_fields[i]->GetPhys(), 1, m_previousSolution[i],
934 1);
935 }
936
937 return false;
938}
939
940/**
941 *
942 */
943void UnsteadySystem::v_SteadyStateResidual([[maybe_unused]] int step,
945{
946 const int nPoints = GetTotPoints();
947 const int nFields = m_fields.size();
948
949 // Holds L2 errors.
950 Array<OneD, NekDouble> RHSL2(nFields);
951 Array<OneD, NekDouble> residual(nFields);
952 Array<OneD, NekDouble> reference(nFields);
953
954 for (int i = 0; i < nFields; ++i)
955 {
956 Array<OneD, NekDouble> tmp(nPoints);
957
958 Vmath::Vsub(nPoints, m_fields[i]->GetPhys(), 1, m_previousSolution[i],
959 1, tmp, 1);
960 Vmath::Vmul(nPoints, tmp, 1, tmp, 1, tmp, 1);
961 residual[i] = Vmath::Vsum(nPoints, tmp, 1);
962
964 tmp, 1);
965 reference[i] = Vmath::Vsum(nPoints, tmp, 1);
966 }
967
968 m_comm->GetSpaceComm()->AllReduce(residual, LibUtilities::ReduceSum);
969 m_comm->GetSpaceComm()->AllReduce(reference, LibUtilities::ReduceSum);
970
971 // L2 error.
972 for (int i = 0; i < nFields; ++i)
973 {
974 reference[i] = (reference[i] == 0) ? 1 : reference[i];
975 L2[i] = sqrt(residual[i] / reference[i]);
976 }
977}
978
979/**
980 *
981 */
983 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
985 [[maybe_unused]] const NekDouble time)
986{
987
988 if (&inarray != &outarray)
989 {
990 for (int i = 0; i < inarray.size(); ++i)
991 {
992 Vmath::Vcopy(GetNpoints(), inarray[i], 1, outarray[i], 1);
993 }
994 }
995}
996
997} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:223
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::string RegisterCmdLineArgument(const std::string &pName, const std::string &pShortName, const std::string &pDescription)
Registers a command-line argument with the session reader.
Binds a set of functions for use by time integration schemes.
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:65
virtual SOLVER_UTILS_EXPORT void v_ALEPreMultiplyMass(Array< OneD, Array< OneD, NekDouble > > &fields)
Definition: ALEHelper.cpp:108
virtual SOLVER_UTILS_EXPORT void v_ALEInitObject(int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Definition: ALEHelper.cpp:42
SOLVER_UTILS_EXPORT void ALEDoElmtInvMass(Array< OneD, Array< OneD, NekDouble > > &traceNormals, Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
Update m_fields with u^n by multiplying by inverse mass matrix. That's then used in e....
Definition: ALEHelper.cpp:131
virtual SOLVER_UTILS_EXPORT void v_UpdateGridVelocity(const NekDouble &time)
Definition: ALEHelper.cpp:90
A base class for describing how to solve specific equations.
int m_spacedim
Spatial dimension (>= expansion dim).
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
int m_infosteps
Number of time steps between outputting status information.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
NekDouble m_fintime
Finish time of the simulation.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
NekDouble m_checktime
Time between checkpoints.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareFeld=true)
Initialisation object for EquationSystem.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
int m_initialStep
Number of the step where the simulation should begin.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT int GetTotPoints()
int m_nchk
Number of checkpoints written so far.
int m_steps
Number of steps to take.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &l)
Virtual function for generating summary information.
int m_checksteps
Number of steps between checkpoints.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_CFLGrowth
CFL growth rate.
bool m_explicitReaction
Indicates if explicit or implicit treatment of reaction is used.
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
Storage for previous solution for steady-state check.
SOLVER_UTILS_EXPORT void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck()
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
SOLVER_UTILS_EXPORT void DoDummyProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Perform dummy projection.
NekDouble m_steadyStateTol
Tolerance to which steady state should be evaluated at.
int m_abortSteps
Number of steps between checks for abort conditions.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator()
Get the maximum timestep estimator for cfl control.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate(int step)
NekDouble m_CFLEnd
Maximun cfl in cfl growth.
virtual SOLVER_UTILS_EXPORT void v_PrintSummaryStatistics(const NekDouble intTime)
Print Summary Statistics.
virtual SOLVER_UTILS_EXPORT void v_PrintStatusInformation(const int step, const NekDouble cpuTime)
Print Status Information.
bool CheckSteadyState(int step, const NekDouble &totCPUTime=0.0)
Calculate whether the system has reached a steady state by observing residuals to a user-defined tole...
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.
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual(int step, Array< OneD, NekDouble > &L2)
SOLVER_UTILS_EXPORT void AppendOutput1D(void)
Print the solution at each solution point in a txt file.
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate(int step)
SOLVER_UTILS_EXPORT void CheckForRestartTime(NekDouble &time, int &nchk)
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans()
SOLVER_UTILS_EXPORT void SteadyStateResidual(int step, Array< OneD, NekDouble > &L2)
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperators & GetTimeIntegrationSchemeOperators()
Returns the time integration scheme operators.
int m_steadyStateSteps
Check for steady state at step interval.
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeSharedPtr & GetTimeIntegrationScheme()
Returns the time integration scheme.
int m_filtersInfosteps
Number of time steps between outputting filters information.
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...
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
SOLVER_UTILS_EXPORT NekDouble GetTimeStep()
SOLVER_UTILS_EXPORT void v_DoSolve() override
Solves an unsteady problem.
TimeIntegrationSchemeFactory & GetTimeIntegrationSchemeFactory()
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:322
static std::string PortablePath(const fs::path &path)
create portable path on different platforms for std::filesystem path.
Definition: Filesystem.hpp:66
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:51
std::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
static const NekDouble kNekZeroTol
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
FilterFactory & GetFilterFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:375
std::vector< double > z(NPUPPER)
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
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.hpp:72
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.hpp:366
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.hpp:743
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.hpp:644
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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.hpp:220
STL namespace.
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285
std::vector< NekDouble > freeParams
Definition: SessionReader.h:93