Nektar++
EquationSystem.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: EquationSystem.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:
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
43
47
50
52
53#include <boost/format.hpp>
54
55#include <iostream>
56#include <string>
57
58using namespace std;
59
60namespace Nektar::SolverUtils
61{
62
65 LibUtilities::SessionReader::RegisterEnumValue("DEALIASING", "False", 1)};
66
68 LibUtilities::SessionReader::RegisterEnumValue("Projection", "Continuous",
70 LibUtilities::SessionReader::RegisterEnumValue("Projection", "CONTINUOUS",
77 "Projection", "DisContinuous", MultiRegions::eDiscontinuous),
79 "Projection", "Mixed_CG_Discontinuous",
82 "Projection", "MixedCGDG", MultiRegions::eMixed_CG_Discontinuous),
83};
84
85/**
86 * @class EquationSystem
87 *
88 * This class is a base class for all solver implementations. It
89 * provides the underlying generic functionality and interface for
90 * solving equations.
91 *
92 * To solve a steady-state equation, create a derived class from this
93 * class and reimplement the virtual functions to provide custom
94 * implementation for the problem.
95 *
96 * To solve unsteady problems, derive from the UnsteadySystem class
97 * instead which provides general time integration.
98 */
100{
101 static EquationSystemFactory instance;
102 return instance;
103}
104
105/**
106 * This constructor is protected as the objects of this class are never
107 * instantiated directly.
108 * @param pSession The session reader holding problem parameters.
109 */
113 : m_comm(pSession->GetComm()), m_session(pSession), m_graph(pGraph),
114 m_lambda(0), m_infosteps(10),
115 m_fieldMetaDataMap(LibUtilities::NullFieldMetaDataMap)
116{
117 // set up session names in fieldMetaDataMap
118 const vector<std::string> filenames = m_session->GetFilenames();
119
120 for (int i = 0; i < filenames.size(); ++i)
121 {
122 string sessionname = "SessionName";
123 sessionname += std::to_string(i);
124 m_fieldMetaDataMap[sessionname] = filenames[i];
125 m_fieldMetaDataMap["ChkFileNum"] = std::to_string(0);
126 }
127}
128
129/**
130 * @brief Destructor for class EquationSystem.
131 */
133{
135 LocalRegions::MatrixKey::opLess>::ClearManager();
137 LocalRegions::MatrixKey::opLess>::ClearManager();
138}
139
140/**
141 * @brief Initialisation object for EquationSystem.
142 */
143void EquationSystem::v_InitObject(bool DeclareFields)
144{
145 // Save the basename of input file name for output details
146 m_sessionName = m_session->GetSessionName();
147
148 // Instantiate a field reader/writer
150
151 // Also read and store the boundary conditions
155
156 // Set space dimension for use in class
157 m_spacedim = m_graph->GetSpaceDimension();
158
159 // Setting parameteres for homogenous problems
160 m_HomoDirec = 0;
161 m_useFFT = false;
162 m_homogen_dealiasing = false;
163 m_singleMode = false;
164 m_halfMode = false;
165 m_multipleModes = false;
167
168 m_verbose = m_session->DefinesCmdLineArgument("verbose");
169
170 if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
171 {
172 std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
173 m_spacedim = 3;
174
175 if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
176 (HomoStr == "1D") || (HomoStr == "Homo1D"))
177 {
179 m_session->LoadParameter("LZ", m_LhomZ);
180 m_HomoDirec = 1;
181
182 if (m_session->DefinesSolverInfo("ModeType"))
183 {
184 m_session->MatchSolverInfo("ModeType", "SingleMode",
185 m_singleMode, false);
186 m_session->MatchSolverInfo("ModeType", "HalfMode", m_halfMode,
187 false);
188 m_session->MatchSolverInfo("ModeType", "MultipleModes",
189 m_multipleModes, false);
190 }
191
192 // Stability Analysis flags
193 if (m_session->DefinesSolverInfo("ModeType"))
194 {
195 if (m_singleMode)
196 {
197 m_npointsZ = 2;
198 }
199 else if (m_halfMode)
200 {
201 m_npointsZ = 1;
202 }
203 else if (m_multipleModes)
204 {
205 m_npointsZ = m_session->GetParameter("HomModesZ");
206 }
207 else
208 {
209 ASSERTL0(false, "SolverInfo ModeType not valid");
210 }
211 }
212 else
213 {
214 m_npointsZ = m_session->GetParameter("HomModesZ");
215 }
216 }
217
218 if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D") ||
219 (HomoStr == "2D") || (HomoStr == "Homo2D"))
220 {
222 m_session->LoadParameter("HomModesY", m_npointsY);
223 m_session->LoadParameter("LY", m_LhomY);
224 m_session->LoadParameter("HomModesZ", m_npointsZ);
225 m_session->LoadParameter("LZ", m_LhomZ);
226 m_HomoDirec = 2;
227 }
228
229 if ((HomoStr == "HOMOGENEOUS3D") || (HomoStr == "Homogeneous3D") ||
230 (HomoStr == "3D") || (HomoStr == "Homo3D"))
231 {
233 m_session->LoadParameter("HomModesY", m_npointsY);
234 m_session->LoadParameter("LY", m_LhomY);
235 m_session->LoadParameter("HomModesZ", m_npointsZ);
236 m_session->LoadParameter("LZ", m_LhomZ);
237 m_HomoDirec = 2;
238 }
239
240 m_session->MatchSolverInfo("USEFFT", "FFTW", m_useFFT, false);
241
242 m_session->MatchSolverInfo("DEALIASING", "True", m_homogen_dealiasing,
243 false);
244 }
245 else
246 {
247 // set to default value so can use to identify 2d or 3D
248 // (homogeneous) expansions
249 m_npointsZ = 1;
250 }
251
252 m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "True",
253 m_specHP_dealiasing, false);
254 if (m_specHP_dealiasing == false)
255 {
256 m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "On",
257 m_specHP_dealiasing, false);
258 }
259
260 // Options to determine type of projection from file or directly
261 // from constructor
262 if (m_session->DefinesSolverInfo("PROJECTION"))
263 {
264 std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
265
266 if ((ProjectStr == "Continuous") || (ProjectStr == "Galerkin") ||
267 (ProjectStr == "CONTINUOUS") || (ProjectStr == "GALERKIN"))
268 {
270 }
271 else if ((ProjectStr == "MixedCGDG") ||
272 (ProjectStr == "Mixed_CG_Discontinuous"))
273 {
275 }
276 else if (ProjectStr == "DisContinuous")
277 {
279 }
280 else
281 {
282 ASSERTL0(false, "PROJECTION value not recognised");
283 }
284 }
285 else
286 {
287 cerr << "Projection type not specified in SOLVERINFO,"
288 "defaulting to continuous Galerkin"
289 << endl;
291 }
292
293 // Enforce singularity check for some problems
295
296 int i;
297 int nvariables = m_session->GetVariables().size();
298 bool DeclareCoeffPhysArrays = true;
299
301 m_spacedim = m_graph->GetSpaceDimension() + m_HomoDirec;
302 m_expdim = m_graph->GetMeshDimension();
303
304 if (DeclareFields) // declare field if required
305 {
306 /// Continuous field
309 {
310 switch (m_expdim)
311 {
312 case 1:
313 {
316 {
317 const LibUtilities::PointsKey PkeyY(
319 const LibUtilities::BasisKey BkeyY(
321 const LibUtilities::PointsKey PkeyZ(
323 const LibUtilities::BasisKey BkeyZ(
325
326 for (i = 0; i < m_fields.size(); i++)
327 {
330 AllocateSharedPtr(m_session, BkeyY, BkeyZ,
333 m_session->GetVariable(i));
334 }
335 }
336 else
337 {
338 for (i = 0; i < m_fields.size(); i++)
339 {
340 m_fields[i] =
344 m_session->GetVariable(i));
345 }
346 }
347 break;
348 }
349 case 2:
350 {
352 {
353 // Fourier single mode stability analysis
354 if (m_singleMode)
355 {
356 const LibUtilities::PointsKey PkeyZ(
359
360 const LibUtilities::BasisKey BkeyZ(
362 PkeyZ);
363
364 for (i = 0; i < m_fields.size(); i++)
365 {
368 AllocateSharedPtr(
369 m_session, BkeyZ, m_LhomZ, m_useFFT,
371 m_session->GetVariable(i),
373 }
374 }
375 // Half mode stability analysis
376 else if (m_halfMode)
377 {
378 const LibUtilities::PointsKey PkeyZ(
381
382 const LibUtilities::BasisKey BkeyZR(
384 PkeyZ);
385
386 const LibUtilities::BasisKey BkeyZI(
388 PkeyZ);
389
390 for (i = 0; i < m_fields.size(); i++)
391 {
392 if (m_session->GetVariable(i).compare("w") == 0)
393 {
395 MultiRegions::
396 ContField3DHomogeneous1D>::
397 AllocateSharedPtr(
398 m_session, BkeyZI, m_LhomZ,
400 m_graph, m_session->GetVariable(i),
402 }
403 else
404 {
406 MultiRegions::
407 ContField3DHomogeneous1D>::
408 AllocateSharedPtr(
409 m_session, BkeyZR, m_LhomZ,
411 m_graph, m_session->GetVariable(i),
413 }
414 }
415 }
416 // Normal homogeneous 1D
417 else
418 {
419 const LibUtilities::PointsKey PkeyZ(
421 const LibUtilities::BasisKey BkeyZ(
423
424 for (i = 0; i < m_fields.size(); i++)
425 {
428 AllocateSharedPtr(
429 m_session, BkeyZ, m_LhomZ, m_useFFT,
431 m_session->GetVariable(i),
433 }
434 }
435 }
436 else
437 {
438 i = 0;
442 m_session->GetVariable(i),
443 DeclareCoeffPhysArrays,
445 m_fields[0] = firstfield;
446 for (i = 1; i < m_fields.size(); i++)
447 {
448 if (m_graph->SameExpansionInfo(
449 m_session->GetVariable(0),
450 m_session->GetVariable(i)))
451 {
452 m_fields[i] =
455 *firstfield, m_graph,
456 m_session->GetVariable(i),
457 DeclareCoeffPhysArrays,
459 }
460 else
461 {
462 m_fields[i] =
466 m_session->GetVariable(i),
467 DeclareCoeffPhysArrays,
469 }
470 }
471
472 if (m_projectionType ==
474 {
475 /// Setting up the normals
478
479 for (i = 0; i < m_spacedim; ++i)
480 {
481 m_traceNormals[i] =
483 }
484
485 m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
486 }
487 }
488
489 break;
490 }
491 case 3:
492 {
493 i = 0;
497 m_session->GetVariable(i),
498 DeclareCoeffPhysArrays,
500
501 m_fields[0] = firstfield;
502 for (i = 1; i < m_fields.size(); i++)
503 {
504 if (m_graph->SameExpansionInfo(
505 m_session->GetVariable(0),
506 m_session->GetVariable(i)))
507 {
508 m_fields[i] =
511 *firstfield, m_graph,
512 m_session->GetVariable(i),
513 DeclareCoeffPhysArrays,
515 }
516 else
517 {
518 m_fields[i] =
522 m_session->GetVariable(i),
523 DeclareCoeffPhysArrays,
525 }
526 }
527
528 if (m_projectionType ==
530 {
531 /// Setting up the normals
534 for (i = 0; i < m_spacedim; ++i)
535 {
536 m_traceNormals[i] =
538 }
539
540 m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
541 // Call the trace on all fields to ensure DG setup.
542 for (i = 1; i < m_fields.size(); ++i)
543 {
544 m_fields[i]->GetTrace();
545 }
546 }
547 break;
548 }
549 default:
550 ASSERTL0(false, "Expansion dimension not recognised");
551 break;
552 }
553 }
554 // Discontinuous field
555 else
556 {
557 switch (m_expdim)
558 {
559 case 1:
560 {
563 {
564 const LibUtilities::PointsKey PkeyY(
566 const LibUtilities::BasisKey BkeyY(
568 const LibUtilities::PointsKey PkeyZ(
570 const LibUtilities::BasisKey BkeyZ(
572
573 for (i = 0; i < m_fields.size(); i++)
574 {
577 AllocateSharedPtr(m_session, BkeyY, BkeyZ,
580 m_session->GetVariable(i));
581 }
582 }
583 else
584 {
585 for (i = 0; i < m_fields.size(); i++)
586 {
587 m_fields[i] =
591 m_session->GetVariable(i));
592 }
593 }
594
595 break;
596 }
597 case 2:
598 {
600 {
601 const LibUtilities::PointsKey PkeyZ(
603 const LibUtilities::BasisKey BkeyZ(
605
606 for (i = 0; i < m_fields.size(); i++)
607 {
610 AllocateSharedPtr(m_session, BkeyZ, m_LhomZ,
611 m_useFFT,
613 m_session->GetVariable(i));
614 }
615 }
616 else
617 {
618 i = 0;
622 m_session->GetVariable(i));
623 m_fields[0] = firstfield;
624 for (i = 1; i < m_fields.size(); i++)
625 {
626 if (m_graph->SameExpansionInfo(
627 m_session->GetVariable(0),
628 m_session->GetVariable(i)))
629 {
630 m_fields[i] =
633 *firstfield, m_graph,
634 m_session->GetVariable(i));
635 }
636 else
637 {
638 m_fields[i] =
642 m_session->GetVariable(i));
643 }
644 }
645 }
646
647 break;
648 }
649 case 3:
650 {
652 {
653 ASSERTL0(
654 false,
655 "3D fully periodic problems not implemented yet");
656 }
657 else
658 {
659 i = 0;
663 m_session->GetVariable(i));
664 m_fields[0] = firstfield;
665 for (i = 1; i < m_fields.size(); i++)
666 {
667 if (m_graph->SameExpansionInfo(
668 m_session->GetVariable(0),
669 m_session->GetVariable(i)))
670 {
671 m_fields[i] =
674 *firstfield, m_graph,
675 m_session->GetVariable(i));
676 }
677 else
678 {
679 m_fields[i] =
683 m_session->GetVariable(i));
684 }
685 }
686 }
687 break;
688 }
689 default:
690 ASSERTL0(false, "Expansion dimension not recognised");
691 break;
692 }
693
694 // Setting up the normals
696
697 for (i = 0; i < m_spacedim; ++i)
698 {
699 m_traceNormals[i] =
701 }
702
703 m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
704 }
705 // Zero all physical fields initially
707 }
708
709 // Set Default Parameter
710 m_session->LoadParameter("Time", m_time, 0.0);
711 m_session->LoadParameter("TimeStep", m_timestep, 0.0);
712 m_session->LoadParameter("NumSteps", m_steps, 0);
713 m_session->LoadParameter("IO_CheckSteps", m_checksteps, 0);
714 m_session->LoadParameter("IO_CheckTime", m_checktime, 0.0);
715 m_session->LoadParameter("FinTime", m_fintime, 0);
716 m_session->LoadParameter("NumQuadPointsError", m_NumQuadPointsError, 0);
717
718 // Check uniqueness of checkpoint output
719 ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
720 (m_checktime > 0.0 && m_checksteps == 0) ||
721 (m_checktime == 0.0 && m_checksteps > 0),
722 "Only one of IO_CheckTime and IO_CheckSteps "
723 "should be set!");
724 m_session->LoadParameter("TimeIncrementFactor", m_TimeIncrementFactor, 1.0);
725
726 // Check for parallel-in-time
727 if (m_comm->IsParallelInTime())
728 {
729 ASSERTL0(m_fintime == 0.0,
730 "Only specify NumSteps and TimeSteps for Parallel-in-Time. "
731 "FinTime should not be used! ");
732 }
733
734 m_nchk = 0;
735 m_iterPIT = 0;
736}
737
738/**
739 *
740 */
742 std::string name, const MultiRegions::ExpListSharedPtr &field, bool cache)
743{
745 if (!field)
746 {
747 vField = m_fields[0];
748 }
749
750 if (cache)
751 {
752 if ((m_sessionFunctions.find(name) == m_sessionFunctions.end()) ||
753 (m_sessionFunctions[name]->GetSession() != m_session) ||
754 (m_sessionFunctions[name]->GetExpansion() != vField))
755 {
758 m_session, vField, name, cache);
759 }
760
761 return m_sessionFunctions[name];
762 }
763 else
764 {
766 new SessionFunction(m_session, vField, name, cache));
767 }
768}
769
770/**
771 * If boundary conditions are time-dependent, they will be evaluated at
772 * the time specified.
773 * @param time The time at which to evaluate the BCs
774 */
776{
777 std::string varName;
778 int nvariables = m_fields.size();
779 for (int i = 0; i < nvariables; ++i)
780 {
781 varName = m_session->GetVariable(i);
782 m_fields[i]->EvaluateBoundaryConditions(time, varName);
783 }
784}
785
786/**
787 * Compute the error in the L2-norm.
788 * @param field The field to compare.
789 * @param exactsoln The exact solution to compare with.
790 * @param Normalised Normalise L2-error.
791 * @returns Error in the L2-norm.
792 */
794 const Array<OneD, NekDouble> &exactsoln,
795 bool Normalised)
796{
797 NekDouble L2error = -1.0;
798
799 if (m_NumQuadPointsError == 0)
800 {
801 if (m_fields[field]->GetPhysState() == false)
802 {
803 m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
804 m_fields[field]->UpdatePhys());
805 }
806
807 // Transform from (phys, wave) -> (phys, phys), for 2.5D
808 bool physwavetrans = false;
809 if (m_fields[field]->GetWaveSpace() == true)
810 {
811 m_fields[field]->HomogeneousBwdTrans(
812 m_fields[field]->GetTotPoints(), m_fields[field]->GetPhys(),
813 m_fields[field]->UpdatePhys());
814 physwavetrans = true;
815 }
816
817 if (exactsoln.size())
818 {
819 L2error =
820 m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
821 }
822 else if (m_session->DefinesFunction("ExactSolution"))
823 {
825
826 GetFunction("ExactSolution")
827 ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
828
829 L2error =
830 m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
831 }
832 else
833 {
834 L2error = m_fields[field]->L2(m_fields[field]->GetPhys());
835 }
836
837 if (Normalised == true)
838 {
840
841 NekDouble Vol = m_fields[field]->Integral(one);
842 L2error = sqrt(L2error * L2error / Vol);
843 }
844
845 // Transform from (phys, phys) -> (phys, wave), for 2.5D
846 if (physwavetrans == true)
847 {
848 m_fields[field]->HomogeneousFwdTrans(
849 m_fields[field]->GetTotPoints(), m_fields[field]->GetPhys(),
850 m_fields[field]->UpdatePhys());
851 physwavetrans = false;
852 }
853 }
854 else
855 {
856 Array<OneD, NekDouble> errNorms(3);
857 errNorms = ErrorExtraPoints(field);
858 L2error = errNorms[0];
859 }
860 return L2error;
861}
862
863/**
864 * Compute the error in the L_inf-norm
865 * @param field The field to compare.
866 * @param exactsoln The exact solution to compare with.
867 * @returns Error in the L_inft-norm.
868 */
870 const Array<OneD, NekDouble> &exactsoln)
871{
872 NekDouble Linferror = -1.0;
873
874 if (m_NumQuadPointsError == 0)
875 {
876 if (m_fields[field]->GetPhysState() == false)
877 {
878 m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
879 m_fields[field]->UpdatePhys());
880 }
881
882 if (exactsoln.size())
883 {
884 Linferror =
885 m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
886 }
887 else if (m_session->DefinesFunction("ExactSolution"))
888 {
890
891 GetFunction("ExactSolution")
892 ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
893
894 Linferror =
895 m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
896 }
897 else
898 {
899 Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys());
900 }
901 }
902 else
903 {
904 Array<OneD, NekDouble> errNorms(3);
905 errNorms = ErrorExtraPoints(field);
906 Linferror = errNorms[1];
907 }
908
909 return Linferror;
910}
911
912/**
913 * Compute the error in the H1-norm.
914 * @param field The field to compare.
915 * @param exactsoln The exact solution to compare with.
916 * @param Normalised Normalise H1-error.
917 * @returns Error in the H1-norm.
918 */
920 const Array<OneD, NekDouble> &exactsoln,
921 bool Normalised)
922{
923 NekDouble error = -1.0;
924
925 if (m_NumQuadPointsError == 0)
926 {
927 if (m_fields[field]->GetPhysState() == false)
928 {
929 m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
930 m_fields[field]->UpdatePhys());
931 }
932
933 if (exactsoln.size())
934 {
935 error = m_fields[field]->H1(m_fields[field]->GetPhys(), exactsoln);
936 }
937 else if (m_session->DefinesFunction("ExactSolution"))
938 {
940
941 GetFunction("ExactSolution")
942 ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
943
944 error = m_fields[field]->H1(m_fields[field]->GetPhys(), exactsoln);
945 }
946 else
947 {
948 error = m_fields[field]->H1(m_fields[field]->GetPhys());
949 }
950
951 if (Normalised == true)
952 {
954
955 NekDouble Vol = m_fields[field]->Integral(one);
956 error = sqrt(error * error / Vol);
957 }
958 }
959 else
960 {
961 Array<OneD, NekDouble> errNorms(3);
962 errNorms = ErrorExtraPoints(field);
963 error = errNorms[2];
964 }
965 return error;
966}
967
968/**
969 * Compute the error in the L2-, Linf- and H1-norm for a larger number of
970 * quadrature points.
971 * @param field The field to compare.
972 * @returns Error in the L2-, Linf- and H1-norm.
973 */
975{
976 int NumModes = GetNumExpModes();
977 Array<OneD, NekDouble> errNorms(3); // Expects 3 norms (L2, Linf, H1)
978
982 LibUtilities::eGaussRadauMAlpha1Beta0);
988 PkeyT1);
990 PkeyT2);
992 PkeyQ1);
994 PkeyQ2);
995
996 LibUtilities::BasisKeyVector Tkeys, Qkeys;
997
998 // make a copy of the ExpansionInfoMap
999 SpatialDomains::ExpansionInfoMap NewExpInfo = m_graph->GetExpansionInfo();
1002 NewExpInfo);
1003
1004 // reset new graph with new keys
1005 Tkeys.push_back(BkeyT1);
1006 Tkeys.push_back(BkeyT2);
1007 m_graph->ResetExpansionInfoToBasisKey(ExpInfo, LibUtilities::eTriangle,
1008 Tkeys);
1009 Qkeys.push_back(BkeyQ1);
1010 Qkeys.push_back(BkeyQ2);
1011 m_graph->ResetExpansionInfoToBasisKey(ExpInfo, LibUtilities::eQuadrilateral,
1012 Qkeys);
1013
1014 MultiRegions::ExpListSharedPtr ErrorExplist =
1016 NewExpInfo);
1017
1018 int ErrorCoordim = ErrorExplist->GetCoordim(0);
1019 int ErrorNq = ErrorExplist->GetTotPoints();
1020
1021 Array<OneD, NekDouble> ErrorXc0(ErrorNq, 0.0);
1022 Array<OneD, NekDouble> ErrorXc1(ErrorNq, 0.0);
1023 Array<OneD, NekDouble> ErrorXc2(ErrorNq, 0.0);
1024
1025 switch (ErrorCoordim)
1026 {
1027 case 1:
1028 ErrorExplist->GetCoords(ErrorXc0);
1029 break;
1030 case 2:
1031 ErrorExplist->GetCoords(ErrorXc0, ErrorXc1);
1032 break;
1033 case 3:
1034 ErrorExplist->GetCoords(ErrorXc0, ErrorXc1, ErrorXc2);
1035 break;
1036 }
1038 m_session->GetFunction("ExactSolution", field);
1039
1040 // Evaluate the exact solution
1041 Array<OneD, NekDouble> ErrorSol(ErrorNq);
1042
1043 exSol->Evaluate(ErrorXc0, ErrorXc1, ErrorXc2, m_time, ErrorSol);
1044
1045 // Calcualte spectral/hp approximation on the quadrature points
1046 // of this new expansion basis
1047 ErrorExplist->BwdTrans(m_fields[field]->GetCoeffs(),
1048 ErrorExplist->UpdatePhys());
1049
1050 errNorms[0] = ErrorExplist->L2(ErrorExplist->GetPhys(), ErrorSol);
1051 errNorms[1] = ErrorExplist->Linf(ErrorExplist->GetPhys(), ErrorSol);
1052 errNorms[2] = ErrorExplist->H1(ErrorExplist->GetPhys(), ErrorSol);
1053
1054 return errNorms;
1055}
1056
1057/**
1058 * Set the physical fields based on a restart file, or a function
1059 * describing the initial condition given in the session.
1060 * @param initialtime Time at which to evaluate the function.
1061 * @param dumpInitialConditions Write the initial condition to file?
1062 */
1064 [[maybe_unused]] NekDouble initialtime, bool dumpInitialConditions,
1065 const int domain)
1066{
1067 if (m_session->GetComm()->GetRank() == 0)
1068 {
1069 cout << "Initial Conditions:" << endl;
1070 }
1071
1072 if (m_session->DefinesFunction("InitialConditions"))
1073 {
1074 GetFunction("InitialConditions")
1075 ->Evaluate(m_session->GetVariables(), m_fields, m_time, domain);
1076 // Enforce C0 Continutiy of initial condiiton
1079 {
1080 for (int i = 0; i < m_fields.size(); ++i)
1081 {
1082 m_fields[i]->LocalToGlobal();
1083 m_fields[i]->GlobalToLocal();
1084 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
1085 m_fields[i]->UpdatePhys());
1086 }
1087 }
1088
1089 if (m_session->GetComm()->GetRank() == 0)
1090 {
1091 for (int i = 0; i < m_fields.size(); ++i)
1092 {
1093 std::string varName = m_session->GetVariable(i);
1094 cout << " - Field " << varName << ": "
1095 << GetFunction("InitialConditions")
1096 ->Describe(varName, domain)
1097 << endl;
1098 }
1099 }
1100 }
1101 else
1102 {
1103 int nq = m_fields[0]->GetNpoints();
1104 for (int i = 0; i < m_fields.size(); i++)
1105 {
1106 Vmath::Zero(nq, m_fields[i]->UpdatePhys(), 1);
1107 m_fields[i]->SetPhysState(true);
1108 Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(),
1109 1);
1110 if (m_session->GetComm()->GetRank() == 0)
1111 {
1112 cout << " - Field " << m_session->GetVariable(i)
1113 << ": 0 (default)" << endl;
1114 }
1115 }
1116 }
1117
1118 if (dumpInitialConditions && m_checksteps && m_nchk == 0 &&
1119 !m_comm->IsParallelInTime())
1120 {
1122 }
1123 else if (dumpInitialConditions && m_nchk == 0 && m_comm->IsParallelInTime())
1124 {
1125 std::string newdir = m_sessionName + ".pit";
1126 if (!fs::is_directory(newdir))
1127 {
1128 fs::create_directory(newdir);
1129 }
1130 if (m_comm->GetTimeComm()->GetRank() == 0)
1131 {
1132 WriteFld(newdir + "/" + m_sessionName + "_0.fld");
1133 }
1134 }
1135 ++m_nchk;
1136}
1137
1138/**
1139 *
1140 */
1142 Array<OneD, NekDouble> &outfield,
1143 const NekDouble time)
1144{
1145 ASSERTL0(outfield.size() == m_fields[field]->GetNpoints(),
1146 "ExactSolution array size mismatch.");
1147 Vmath::Zero(outfield.size(), outfield, 1);
1148 if (m_session->DefinesFunction("ExactSolution"))
1149 {
1150 GetFunction("ExactSolution")
1151 ->Evaluate(m_session->GetVariable(field), outfield, time);
1152 }
1153}
1154
1155/**
1156 * By default, nothing needs initialising at the EquationSystem level.
1157 */
1158void EquationSystem::v_DoInitialise([[maybe_unused]] bool dumpInitialConditions)
1159{
1160}
1161
1162/**
1163 *
1164 */
1166{
1167}
1168
1169/**
1170 * Virtual function to define if operator in DoSolve is
1171 * negated with regard to the strong form. This is currently
1172 * only used in Arnoldi solves. Default is false.
1173 */
1175{
1176 return false;
1177}
1178
1179/**
1180 *
1181 */
1183{
1184}
1185
1186/**
1187 *
1188 */
1190{
1191}
1192
1193/// Virtual function for generating summary information.
1195{
1196 SessionSummary(l);
1197}
1198
1199/**
1200 * Write the field data to file. The file is named according to the
1201 * session name with the extension .fld appended.
1202 */
1204{
1205 if (!m_comm->IsParallelInTime())
1206 {
1207 // Serial-in-time
1208 WriteFld(m_sessionName + ".fld");
1209 }
1210 else
1211 {
1212 // Parallel-in-time
1213 std::string newdir = m_sessionName + ".pit";
1214 if (!fs::is_directory(newdir))
1215 {
1216 fs::create_directory(newdir);
1217 }
1218 WriteFld(newdir + "/" + m_sessionName + "_" +
1219 std::to_string(m_windowPIT * m_comm->GetTimeComm()->GetSize() +
1220 m_comm->GetTimeComm()->GetRank() + 1) +
1221 ".fld");
1222 }
1223}
1224
1225/**
1226 * Zero the physical fields.
1227 */
1229{
1230 for (int i = 0; i < m_fields.size(); i++)
1231 {
1232 Vmath::Zero(m_fields[i]->GetNpoints(), m_fields[i]->UpdatePhys(), 1);
1233 }
1234}
1235
1236/**
1237 * FwdTrans the m_fields members
1238 */
1240{
1241 for (int i = 0; i < m_fields.size(); i++)
1242 {
1243 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1244 m_fields[i]->UpdateCoeffs());
1245 m_fields[i]->SetPhysState(false);
1246 }
1247}
1248
1249/**
1250 * Write the n-th checkpoint file.
1251 * @param n The index of the checkpoint file.
1252 */
1254{
1255 if (!m_comm->IsParallelInTime())
1256 {
1257 // Serial-in-time
1258 std::string outname = m_sessionName + "_" + std::to_string(n);
1259 WriteFld(outname + ".chk");
1260 }
1261 else
1262 {
1263 // Parallel-in-time
1264 auto loc = m_sessionName.find("_xml/");
1265 auto sessionName = m_sessionName.substr(0, loc);
1266 std::string paradir =
1267 sessionName + "_" + std::to_string(m_iterPIT) + ".pit";
1268 if (!fs::is_directory(paradir))
1269 {
1270 fs::create_directory(paradir);
1271 }
1272 std::string outname = paradir + "/" + sessionName + "_timeLevel" +
1273 std::to_string(m_session->GetTimeLevel()) + "_" +
1274 std::to_string(n);
1275 WriteFld(outname + ".chk");
1276 }
1277}
1278
1279/**
1280 * Write the n-th checkpoint file.
1281 * @param n The index of the checkpoint file.
1282 */
1285 std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
1286 std::vector<std::string> &variables)
1287{
1288 if (!m_comm->IsParallelInTime())
1289 {
1290 // Serial-in-time
1291 std::string outname = m_sessionName + "_" + std::to_string(n);
1292 WriteFld(outname, field, fieldcoeffs, variables);
1293 }
1294 else
1295 {
1296 ASSERTL0(false, "Not Implemented for Parallel-in-Time");
1297 }
1298}
1299
1300/**
1301 * Write the n-th base flow into a .chk file
1302 * @param n The index of the base flow file.
1303 */
1305{
1306 std::string outname = m_sessionName + "_BaseFlow_" + std::to_string(n);
1307
1308 WriteFld(outname + ".chk");
1309}
1310
1311/**
1312 * Writes the field data to a file with the given filename.
1313 * @param outname Filename to write to.
1314 */
1315void EquationSystem::WriteFld(const std::string &outname)
1316{
1317 std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_fields.size());
1318 std::vector<std::string> variables(m_fields.size());
1319
1320 for (int i = 0; i < m_fields.size(); ++i)
1321 {
1322 if (m_fields[i]->GetNcoeffs() == m_fields[0]->GetNcoeffs())
1323 {
1324 fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
1325 }
1326 else
1327 {
1328 fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
1329 m_fields[0]->ExtractCoeffsToCoeffs(
1330 m_fields[i], m_fields[i]->GetCoeffs(), fieldcoeffs[i]);
1331 }
1332 variables[i] = m_boundaryConditions->GetVariable(i);
1333 }
1334
1335 ExtraFldOutput(fieldcoeffs, variables);
1336
1337 WriteFld(outname, m_fields[0], fieldcoeffs, variables);
1338}
1339
1340/**
1341 * Writes the field data to a file with the given filename.
1342 * @param outname Filename to write to.
1343 * @param field ExpList on which data is based.
1344 * @param fieldcoeffs An array of array of expansion coefficients.
1345 * @param variables An array of variable names.
1346 */
1347void EquationSystem::WriteFld(const std::string &outname,
1349 std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
1350 std::vector<std::string> &variables)
1351{
1352 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
1353 field->GetFieldDefinitions();
1354 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
1355
1356 // Copy Data into FieldData and set variable
1357 for (int j = 0; j < fieldcoeffs.size(); ++j)
1358 {
1359 for (int i = 0; i < FieldDef.size(); ++i)
1360 {
1361 // Could do a search here to find correct variable
1362 FieldDef[i]->m_fields.push_back(variables[j]);
1363 field->AppendFieldData(FieldDef[i], FieldData[i], fieldcoeffs[j]);
1364 }
1365 }
1366
1367 // Update time in field info if required
1368 if (m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
1369 {
1370 m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
1371 }
1372
1373 // Update step in field info if required
1374 if (m_fieldMetaDataMap.find("ChkFileNum") != m_fieldMetaDataMap.end())
1375 {
1376 m_fieldMetaDataMap["ChkFileNum"] = std::to_string(m_nchk);
1377 }
1378
1379 // If necessary, add mapping information to metadata
1380 // and output mapping coordinates
1382 fields[0] = field;
1386 mapping->Output(fieldMetaDataMap, outname);
1387
1388 // If necessary, add informaton for moving frame reference to metadata
1389 // X, Y, Z translational displacements
1390 // Theta_x, Theta_y, Theta_z angular displacements
1391 // U, V, W translational velocity
1392 // Omega_x, Omega_y, Omega_z angular velocity
1393 // A_x, A_y, A_z translational acceleration
1394 // DOmega_x, DOmega_y, DOmega_z angular acceleration
1395 // X0, Y0, Z0 pivot point
1396 for (size_t i = 0;
1397 i < m_strFrameData.size() && i < m_movingFrameData.size(); ++i)
1398 {
1399 fieldMetaDataMap[m_strFrameData[i]] =
1400 boost::lexical_cast<std::string>(m_movingFrameData[i]);
1401 }
1402
1403 m_fld->Write(outname, FieldDef, FieldData, fieldMetaDataMap,
1404 m_session->GetBackups());
1405}
1406
1407/**
1408 * Import field from infile and load into \a m_fields. This routine will
1409 * also perform a \a BwdTrans to ensure data is in both the physical and
1410 * coefficient storage.
1411 * @param infile Filename to read.
1412 */
1414 const std::string &infile,
1416{
1417 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1418 std::vector<std::vector<NekDouble>> FieldData;
1421 field_fld->Import(infile, FieldDef, FieldData);
1422
1423 // Copy FieldData into m_fields
1424 for (int j = 0; j < pFields.size(); ++j)
1425 {
1426 Vmath::Zero(pFields[j]->GetNcoeffs(), pFields[j]->UpdateCoeffs(), 1);
1427
1428 for (int i = 0; i < FieldDef.size(); ++i)
1429 {
1430 ASSERTL1(FieldDef[i]->m_fields[j] == m_session->GetVariable(j),
1431 std::string("Order of ") + infile +
1432 std::string(" data and that defined in "
1433 "m_boundaryconditions differs"));
1434
1435 pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1436 FieldDef[i]->m_fields[j],
1437 pFields[j]->UpdateCoeffs());
1438 }
1439 pFields[j]->BwdTrans(pFields[j]->GetCoeffs(), pFields[j]->UpdatePhys());
1440 }
1441}
1442
1443/**
1444 * Import field from infile and load into \a m_fields. This routine will
1445 * also perform a \a BwdTrans to ensure data is in both the physical and
1446 * coefficient storage.
1447 * @param infile Filename to read.
1448 * If optionan \a ndomains is specified it assumes we loop over nodmains
1449 * for each nvariables.
1450 */
1452 const std::string &infile,
1453 Array<OneD, MultiRegions::ExpListSharedPtr> &pFields, const int ndomains)
1454{
1455 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1456 std::vector<std::vector<NekDouble>> FieldData;
1457
1458 LibUtilities::Import(infile, FieldDef, FieldData);
1459
1460 int nvariables = GetNvariables();
1461
1462 ASSERTL0(
1463 ndomains * nvariables == pFields.size(),
1464 "Number of fields does not match the number of variables and domains");
1465
1466 // Copy FieldData into m_fields
1467 for (int j = 0; j < ndomains; ++j)
1468 {
1469 for (int i = 0; i < nvariables; ++i)
1470 {
1471 Vmath::Zero(pFields[j * nvariables + i]->GetNcoeffs(),
1472 pFields[j * nvariables + i]->UpdateCoeffs(), 1);
1473
1474 for (int n = 0; n < FieldDef.size(); ++n)
1475 {
1476 ASSERTL1(FieldDef[n]->m_fields[i] == m_session->GetVariable(i),
1477 std::string("Order of ") + infile +
1478 std::string(" data and that defined in "
1479 "m_boundaryconditions differs"));
1480
1481 pFields[j * nvariables + i]->ExtractDataToCoeffs(
1482 FieldDef[n], FieldData[n], FieldDef[n]->m_fields[i],
1483 pFields[j * nvariables + i]->UpdateCoeffs());
1484 }
1485 pFields[j * nvariables + i]->BwdTrans(
1486 pFields[j * nvariables + i]->GetCoeffs(),
1487 pFields[j * nvariables + i]->UpdatePhys());
1488 }
1489 }
1490}
1491
1492/**
1493 * Import field from infile and load into \a pField. This routine will
1494 * also perform a \a BwdTrans to ensure data is in both the physical and
1495 * coefficient storage.
1496 */
1497void EquationSystem::ImportFld(const std::string &infile,
1499 std::string &pFieldName)
1500{
1501 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1502 std::vector<std::vector<NekDouble>> FieldData;
1503
1506 field_fld->Import(infile, FieldDef, FieldData);
1507 int idx = -1;
1508
1509 Vmath::Zero(pField->GetNcoeffs(), pField->UpdateCoeffs(), 1);
1510
1511 for (int i = 0; i < FieldDef.size(); ++i)
1512 {
1513 // find the index of the required field in the file.
1514 for (int j = 0; j < FieldData.size(); ++j)
1515 {
1516 if (FieldDef[i]->m_fields[j] == pFieldName)
1517 {
1518 idx = j;
1519 }
1520 }
1521 ASSERTL1(idx >= 0, "Field " + pFieldName + " not found.");
1522
1523 pField->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1524 FieldDef[i]->m_fields[idx],
1525 pField->UpdateCoeffs());
1526 }
1527 pField->BwdTrans(pField->GetCoeffs(), pField->UpdatePhys());
1528}
1529
1530/**
1531 * Import field from infile and load into the array \a coeffs.
1532 *
1533 * @param infile Filename to read.
1534 * @param fieldStr an array of string identifying fields to be imported
1535 * @param coeffs array of array of coefficients to store imported data
1536 */
1537void EquationSystem::ImportFld(const std::string &infile,
1538 std::vector<std::string> &fieldStr,
1540{
1541
1542 ASSERTL0(fieldStr.size() <= coeffs.size(),
1543 "length of fieldstr should be the same as coeffs");
1544
1545 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1546 std::vector<std::vector<NekDouble>> FieldData;
1547
1550 field_fld->Import(infile, FieldDef, FieldData);
1551
1552 // Copy FieldData into m_fields
1553 for (int j = 0; j < fieldStr.size(); ++j)
1554 {
1555 Vmath::Zero(coeffs[j].size(), coeffs[j], 1);
1556 for (int i = 0; i < FieldDef.size(); ++i)
1557 {
1558 m_fields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1559 fieldStr[j], coeffs[j]);
1560 }
1561 }
1562}
1563
1564/**
1565 * Write out a summary of the session data.
1566 * @param out Output stream to write data to.
1567 */
1569{
1570 if (m_session->DefinesSolverInfo("EQTYPE"))
1571 {
1572 AddSummaryItem(s, "EquationType", m_session->GetSolverInfo("EQTYPE"));
1573 }
1574 AddSummaryItem(s, "Session Name", m_sessionName);
1575 AddSummaryItem(s, "Spatial Dim.", m_spacedim);
1576 AddSummaryItem(s, "Max SEM Exp. Order",
1577 m_fields[0]->EvalBasisNumModesMax());
1578
1579 if (m_session->GetComm()->GetSize() > 1)
1580 {
1581 AddSummaryItem(s, "Num. Processes", m_session->GetComm()->GetSize());
1582 }
1583
1585 {
1586 AddSummaryItem(s, "Quasi-3D", "Homogeneous in z-direction");
1587 AddSummaryItem(s, "Expansion Dim.", m_expdim + 1);
1588 AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
1589 AddSummaryItem(s, "Hom. length (LZ)", m_LhomZ);
1590 AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
1591 if (m_halfMode)
1592 {
1593 AddSummaryItem(s, "ModeType", "Half Mode");
1594 }
1595 else if (m_singleMode)
1596 {
1597 AddSummaryItem(s, "ModeType", "Single Mode");
1598 }
1599 else if (m_multipleModes)
1600 {
1601 AddSummaryItem(s, "ModeType", "Multiple Modes");
1602 }
1603 }
1605 {
1606 AddSummaryItem(s, "Quasi-3D", "Homogeneous in yz-plane");
1607 AddSummaryItem(s, "Expansion Dim.", m_expdim + 2);
1608 AddSummaryItem(s, "Num. Hom. Modes (y)", m_npointsY);
1609 AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
1610 AddSummaryItem(s, "Hom. length (LY)", m_LhomY);
1611 AddSummaryItem(s, "Hom. length (LZ)", m_LhomZ);
1612 AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
1613 }
1614 else
1615 {
1616 AddSummaryItem(s, "Expansion Dim.", m_expdim);
1617 }
1618
1619 if (m_session->DefinesSolverInfo("UpwindType"))
1620 {
1621 AddSummaryItem(s, "Riemann Solver",
1622 m_session->GetSolverInfo("UpwindType"));
1623 }
1624
1625 if (m_session->DefinesSolverInfo("AdvectionType"))
1626 {
1627 std::string AdvectionType;
1628 AdvectionType = m_session->GetSolverInfo("AdvectionType");
1630 s, "Advection Type",
1631 GetAdvectionFactory().GetClassDescription(AdvectionType));
1632 }
1633
1634 if (m_session->DefinesSolverInfo("DiffusionType"))
1635 {
1636 std::string DiffusionType;
1637 DiffusionType = m_session->GetSolverInfo("DiffusionType");
1639 s, "Diffusion Type",
1640 GetDiffusionFactory().GetClassDescription(DiffusionType));
1641 }
1642
1644 {
1645 AddSummaryItem(s, "Projection Type", "Continuous Galerkin");
1646 }
1648 {
1649 AddSummaryItem(s, "Projection Type", "Discontinuous Galerkin");
1650 }
1652 {
1653 AddSummaryItem(s, "Projection Type",
1654 "Mixed Continuous Galerkin and Discontinuous");
1655 }
1656}
1657
1658/**
1659 *
1660 */
1662{
1663 return Array<OneD, bool>(m_session->GetVariables().size(), false);
1664}
1665
1666/**
1667 *
1668 */
1670{
1671 ASSERTL0(false, "This function is not valid for the Base class");
1673 return null;
1674}
1675
1676/**
1677 *
1678 */
1680 [[maybe_unused]] std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
1681 [[maybe_unused]] std::vector<std::string> &variables)
1682{
1683}
1684
1685} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:264
Describes the specification for a Basis.
Definition: Basis.h:45
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
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:194
Provides a generic Factory class.
Defines a specification for a set of points.
Definition: Points.h:50
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
int m_spacedim
Spatial dimension (>= expansion dim).
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
int m_expdim
Expansion dimension.
LibUtilities::FieldIOSharedPtr m_fld
Field input/output.
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
virtual SOLVER_UTILS_EXPORT ~EquationSystem()
Destructor.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
SOLVER_UTILS_EXPORT void ImportFld(const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Input field data from the given file.
NekDouble m_time
Current time of simulation.
int m_iterPIT
Number of parallel-in-time time iteration.
bool m_multipleModes
Flag to determine if use multiple homogenenous modes are used.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
int m_windowPIT
Index of windows for parallel-in-time time iteration.
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff()
Virtual function for transformation to coefficient space.
SOLVER_UTILS_EXPORT int GetNpoints()
std::vector< std::string > m_strFrameData
variable name in m_movingFrameData
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp(void)
Virtual function to identify if operator is negated in DoSolve.
NekDouble m_fintime
Finish time of the simulation.
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
Virtual function for the L_2 error computation between fields and a given exact solution.
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure(void)
static std::string projectionTypeLookupIds[]
SOLVER_UTILS_EXPORT int GetNcoeffs()
int m_npointsZ
number of points in Z direction (if homogeneous)
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
Map of known SessionFunctions.
NekDouble m_checktime
Time between checkpoints.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareFeld=true)
Initialisation object for EquationSystem.
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Virtual function for the L_inf error computation between fields and a given exact solution.
static std::string equationSystemTypeLookupIds[]
SOLVER_UTILS_EXPORT int GetNumExpModes()
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
virtual SOLVER_UTILS_EXPORT void v_DoInitialise(bool dumpInitialConditions=true)
Virtual function for initialisation implementation.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains(const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
Input field data from the given file to multiple domains.
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
int m_npointsY
number of points in Y direction (if homogeneous)
std::string m_sessionName
Name of the session.
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow(const int n)
Write base flow file of m_fields.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
bool m_singleMode
Flag to determine if single homogeneous mode is used.
SOLVER_UTILS_EXPORT void ZeroPhysFields()
int m_HomoDirec
number of homogenous directions
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT void ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints(unsigned int field)
Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf].
Array< OneD, bool > m_checkIfSystemSingular
Flag to indicate if the fields should be checked for singularity.
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
SOLVER_UTILS_EXPORT int GetTotPoints()
int m_nchk
Number of checkpoints written so far.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
virtual SOLVER_UTILS_EXPORT void v_DoSolve()
Virtual function for solve implementation.
virtual SOLVER_UTILS_EXPORT Array< OneD, bool > v_GetSystemSingularChecks()
int m_steps
Number of steps to take.
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
SOLVER_UTILS_EXPORT EquationSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises EquationSystem class members.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Array< OneD, NekDouble > m_movingFrameData
Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z,...
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &l)
Virtual function for generating summary information.
bool m_halfMode
Flag to determine if half homogeneous mode is used.
SOLVER_UTILS_EXPORT void FwdTransFields()
virtual SOLVER_UTILS_EXPORT NekDouble v_H1Error(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
Virtual function for the H_1 error computation between fields and a given exact solution.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
virtual SOLVER_UTILS_EXPORT void v_Output(void)
SOLVER_UTILS_EXPORT int GetNvariables()
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys()
Virtual function for transformation to physical space.
GLOBAL_MAPPING_EXPORT typedef std::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:57
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:287
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:322
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:50
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:51
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:74
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eFourierSingleModeSpaced
1D Non Evenly-spaced points for Single Mode analysis
Definition: PointsType.h:75
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
@ eFourierSingleMode
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:64
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
@ eFourierHalfModeIm
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:68
@ eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:66
@ eFourier
Fourier Expansion .
Definition: BasisType.h:55
std::shared_ptr< DisContField > DisContFieldSharedPtr
Definition: DisContField.h:346
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:39
EquationSystemFactory & GetEquationSystemFactory()
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< SessionFunction > SessionFunctionSharedPtr
std::shared_ptr< ExpansionInfoMap > ExpansionInfoMapShPtr
Definition: MeshGraph.h:143
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:141
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
Definition: NekTypeDefs.hpp:68
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
double NekDouble
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
STL namespace.
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285
Used to lookup the create function in NekManager.
Definition: MatrixKey.h:69