31#if defined(MAIA_GCC_COMPILER)
32#pragma GCC diagnostic push
50template <MInt nDim,
class SysEqn>
54 m_geometry(geometry_),
55 m_sponge(solverId_, comm),
80 if(gridProxy_.isActive()) {
97template <MInt nDim,
class SysEqn>
100 RECORD_TIMER_START(m_timers[Timers::Destructor]);
101 RECORD_TIMER_STOP(m_timers[Timers::Destructor]);
102 RECORD_TIMER_STOP(m_timers[Timers::SolverType]);
112template <MInt nDim,
class SysEqn>
117 std::fill(m_timers.begin(), m_timers.end(), -1);
120 NEW_TIMER_GROUP_NOCREATE(m_timerGroup,
"DgCartesianSolver (solverId = " + to_string(m_solverId) +
")");
121 NEW_TIMER_NOCREATE(m_timers[Timers::SolverType],
"total object lifetime", m_timerGroup);
122 RECORD_TIMER_START(m_timers[Timers::SolverType]);
125 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Constructor],
"Constructor", m_timers[Timers::SolverType]);
126 RECORD_TIMER_START(m_timers[Timers::Constructor]);
129 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Run],
"run", m_timers[Timers::SolverType]);
130 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::RunInit],
"run initialization", m_timers[Timers::Run]);
131 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::InitSolverObjects],
"initialize solver infrastructure",
132 m_timers[Timers::RunInit]);
133 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::InitMortarProjection],
"initMortarProjection",
134 m_timers[Timers::InitSolverObjects]);
135 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::InitData],
"initialize solver data", m_timers[Timers::RunInit]);
136 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::InitialCondition],
"initialCondition", m_timers[Timers::InitData]);
137 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::InitMainLoop],
"initialize main loop", m_timers[Timers::RunInit]);
138 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::MainLoop],
"main loop", m_timers[Timers::Run]);
139 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcTimeStep],
"calcTimeStep", m_timers[Timers::MainLoop]);
140 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::ResetExternalSources],
"resetExternalSources", m_timers[Timers::MainLoop]);
141 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::RungeKuttaStep],
"timeStepRk", m_timers[Timers::MainLoop]);
142 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::TimeDeriv],
"calcDgTimeDerivative", m_timers[Timers::RungeKuttaStep]);
143 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::ResetRHS],
"resetRHS", m_timers[Timers::TimeDeriv]);
144 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Prolong],
"prolong to surfaces", m_timers[Timers::TimeDeriv]);
145 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::ForwardProjection],
"forward projection", m_timers[Timers::TimeDeriv]);
146 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SurfExchange],
"MPI surface exchange", m_timers[Timers::TimeDeriv]);
147 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SurfExchangeComm],
"communication", m_timers[Timers::SurfExchange]);
148 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SECommSend],
"send", m_timers[Timers::SurfExchangeComm]);
149 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SECommRecv],
"recv", m_timers[Timers::SurfExchangeComm]);
150 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SurfExchangeCopy],
"copy operations", m_timers[Timers::SurfExchange]);
151 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SECopySend],
"send", m_timers[Timers::SurfExchangeCopy]);
152 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SECopyRecv],
"recv", m_timers[Timers::SurfExchangeCopy]);
153 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SurfExchangeWait],
"waiting", m_timers[Timers::SurfExchange]);
154 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SEWaitSend],
"send", m_timers[Timers::SurfExchangeWait]);
155 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SEWaitRecv],
"recv", m_timers[Timers::SurfExchangeWait]);
156 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::VolInt],
"calcVolumeIntegral", m_timers[Timers::TimeDeriv]);
157 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Flux],
"flux calculation", m_timers[Timers::TimeDeriv]);
158 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::FluxBndry],
"calcBoundarySurfaceFlux", m_timers[Timers::Flux]);
159 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::FluxInner],
"calcInnerSurfaceFlux", m_timers[Timers::Flux]);
160 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::FluxMPI],
"calcMpiSurfaceFlux", m_timers[Timers::Flux]);
161 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SurfInt],
"surface integrals", m_timers[Timers::TimeDeriv]);
162 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Jacobian],
"applyJacobian", m_timers[Timers::TimeDeriv]);
163 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Sources],
"calcSourceTerms", m_timers[Timers::TimeDeriv]);
164 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::ExternalSources],
"applyExternalSourceTerms", m_timers[Timers::TimeDeriv]);
165 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Sponge],
"calcSpongeTerms", m_timers[Timers::TimeDeriv]);
166 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::TimeInt],
"time integration", m_timers[Timers::RungeKuttaStep]);
167 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::MainLoopIO],
"I/O", m_timers[Timers::MainLoop]);
168 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Analysis],
"solution analysis", m_timers[Timers::MainLoop]);
169 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CleanUp],
"cleanUp", m_timers[Timers::Run]);
170 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Destructor],
"Destructor", m_timers[Timers::SolverType]);
171 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::AdaptiveRefinement],
"adaptive refinement", m_timers[Timers::MainLoop]);
174 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Accumulated],
"selected accumulated timers", m_timers[Timers::SolverType]);
175 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::IO],
"IO", m_timers[Timers::Accumulated]);
176 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SaveSolutionFile],
"saveSolutionFile", m_timers[Timers::IO]);
177 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SaveRestartFile],
"saveRestartFile", m_timers[Timers::IO]);
178 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::AnalyzeTimeStep],
"analyzeTimeStep", m_timers[Timers::Accumulated]);
179 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::MPI],
"MPI", m_timers[Timers::Accumulated]);
180 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::MPIComm],
"communication", m_timers[Timers::MPI]);
181 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::MPICopy],
"copy operations", m_timers[Timers::MPI]);
182 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::MPIWait],
"waiting", m_timers[Timers::MPI]);
185 m_isInitTimers =
true;
196template <MInt nDim,
class SysEqn>
215 m_maxNoSurfaces = -1;
216 m_maxNoSurfaces = Context::getSolverProperty<MInt>(
"maxNoSurfaces", m_solverId, AT_, &m_maxNoSurfaces);
231 m_timeSteps = Context::getSolverProperty<MInt>(
"timeSteps", m_solverId, AT_);
242template <MInt nDim,
class SysEqn>
260 m_alwaysSaveFinalSolution = m_solutionInterval > 0;
261 m_alwaysSaveFinalSolution =
262 Context::getSolverProperty<MBool>(
"alwaysSaveFinalSolution", m_solverId, AT_, &m_alwaysSaveFinalSolution);
278 m_alwaysSaveFinalRestart = m_restartInterval > 0;
279 m_alwaysSaveFinalRestart =
280 Context::getSolverProperty<MBool>(
"alwaysSaveFinalRestart", m_solverId, AT_, &m_alwaysSaveFinalRestart);
296 m_saveNodeVariablesToSolutionFile =
false;
297 m_saveNodeVariablesToSolutionFile = Context::getSolverProperty<MBool>(
"saveNodeVariablesToSolutionFile", m_solverId,
298 AT_, &m_saveNodeVariablesToSolutionFile);
314 m_analysisInterval = 0;
316 m_analysisInterval = Context::getSolverProperty<MInt>(
"analysisInterval", m_solverId, AT_, &m_analysisInterval);
317 if(m_analysisInterval < 0) {
318 TERMM(1,
"Analysis interval must be >= 0 (is: " + to_string(m_analysisInterval) +
").");
322 m_pointData.setInputOutputProperties();
324 m_surfaceData.setInputOutputProperties();
326 m_volumeData.setInputOutputProperties();
329 m_slice.setProperties();
360 m_writeSpongeEta =
false;
361 m_writeSpongeEta = Context::getSolverProperty<MBool>(
"writeSpongeEta", m_solverId, AT_, &m_writeSpongeEta);
372 m_writeTimeDerivative = 0;
373 m_writeTimeDerivative =
374 Context::getSolverProperty<MBool>(
"writeTimeDerivative", m_solverId, AT_, &m_writeTimeDerivative);
386 m_noMinutesEndAutoSave = 10;
387 m_noMinutesEndAutoSave =
388 Context::getSolverProperty<MInt>(
"noMinutesEndAutoSave", m_solverId, AT_, &m_noMinutesEndAutoSave);
398 m_endAutoSaveCheckInterval = 50;
399 m_endAutoSaveCheckInterval =
400 Context::getSolverProperty<MInt>(
"endAutoSaveCheckInterval", m_solverId, AT_, &m_endAutoSaveCheckInterval);
432 m_weightDgRbcElements =
false;
433 m_weightDgRbcElements =
434 Context::getSolverProperty<MBool>(
"weightDgRbcElements", m_solverId, AT_, &m_weightDgRbcElements);
450 m_weightPerNode = 1.0;
451 m_weightPerNode = Context::getSolverProperty<MFloat>(
"weightPerNode", m_solverId, AT_, &m_weightPerNode);
465 m_weightPerElement = 0.0;
466 m_weightPerElement = Context::getSolverProperty<MFloat>(
"weightPerElement", m_solverId, AT_, &m_weightPerElement);
498 m_writeInitialSolution =
true;
499 m_writeInitialSolution =
500 Context::getSolverProperty<MBool>(
"writeInitialSolution", m_solverId, AT_, &m_writeInitialSolution);
514 m_writeNodeVarsFile =
true;
515 m_writeNodeVarsFile = Context::getSolverProperty<MBool>(
"writeNodeVarsFile", m_solverId, AT_, &m_writeNodeVarsFile);
525template <MInt nDim,
class SysEqn>
541 m_startTime = Context::getSolverProperty<MFloat>(
"startTime", m_solverId, AT_);
555 m_finalTime = Context::getSolverProperty<MFloat>(
"finalTime", m_solverId, AT_);
557 if(m_startTime > m_finalTime || m_startTime < F0) {
558 cerr <<
"Illegal setup of the time integration parameters. Aborting..." << endl;
559 m_log <<
"Illegal setup of the time integration parameters. Aborting..." << endl;
560 mTerm(1, AT_,
"Illegal setup of the time integration parameters.");
577 m_calcTimeStepInterval = 0;
578 m_calcTimeStepInterval =
579 Context::getSolverProperty<MInt>(
"calcTimeStepInterval", m_solverId, AT_, &m_calcTimeStepInterval);
580 if(m_calcTimeStepInterval < 0) {
581 TERMM(1,
"Recalculation interval for time step must be >= 0 (is: " + to_string(m_calcTimeStepInterval) +
").");
593 m_sbpMode = Context::getSolverProperty<MBool>(
"sbpMode", m_solverId, AT_, &m_sbpMode);
604 m_sbpOperator = Context::getSolverProperty<MString>(
"sbpOperator", m_solverId, AT_, &m_sbpOperator);
620 MString dgIntegrationMethod = m_sbpMode ?
"DG_INTEGRATE_GAUSS_LOBATTO" :
"DG_INTEGRATE_GAUSS";
621 m_dgIntegrationMethod =
622 string2enum(Context::getSolverProperty<MString>(
"dgIntegrationMethod", m_solverId, AT_, &dgIntegrationMethod));
641 MString dgTimeIntegrationScheme =
"DG_TIMEINTEGRATION_CARPENTER_4_5";
643 Context::getSolverProperty<MString>(
"dgTimeIntegrationScheme", m_solverId, AT_, &dgTimeIntegrationScheme));
647 switch(m_dgTimeIntegrationScheme) {
652 m_timeIntegrationCoefficientsA.resize(m_noRkStages);
653 m_timeIntegrationCoefficientsB.resize(m_noRkStages);
654 m_timeIntegrationCoefficientsC.resize(m_noRkStages);
656 m_timeIntegrationCoefficientsA = {{F0, 567301805773.0 / 1357537059087.0, 2404267990393.0 / 2016746695238.0,
657 3550918686646.0 / 2091501179385.0, 1275806237668.0 / 842570457699.0}};
658 m_timeIntegrationCoefficientsB = {{1432997174477.0 / 9575080441755.0, 5161836677717.0 / 13612068292357.0,
659 1720146321549.0 / 2090206949498.0, 3134564353537.0 / 4481467310338.0,
660 2277821191437.0 / 14882151754819.0}};
661 m_timeIntegrationCoefficientsC = {{F0, 1432997174477.0 / 9575080441755.0, 2526269341429.0 / 6820363962896.0,
662 2006345519317.0 / 3224310063776.0, 2802321613138.0 / 2924317926251.0}};
669 m_timeIntegrationCoefficientsA.resize(m_noRkStages);
670 m_timeIntegrationCoefficientsB.resize(m_noRkStages);
671 m_timeIntegrationCoefficientsC.resize(m_noRkStages);
673 m_timeIntegrationCoefficientsA = {{F0, 0.7212962482279240, 0.0107733657161298, 0.5162584698930970,
674 1.7301002866322010, 5.2001293044030760, -0.7837058945416420,
675 0.5445836094332190}};
676 m_timeIntegrationCoefficientsB = {{0.2165936736758085, 0.1773950826411583, 0.0180253861162329, 0.0847347637254149,
677 0.8129106974622483, 1.9034160304227600, 0.1314841743399048,
678 0.2082583170674149}};
679 m_timeIntegrationCoefficientsC = {{F0, 0.2165936736758085, 0.2660343487538170, 0.2840056122522720,
680 0.3251266843788570, 0.4555149599187530, 0.7713219317101170,
681 0.9199028964538660}};
688 m_timeIntegrationCoefficientsA.resize(m_noRkStages);
689 m_timeIntegrationCoefficientsB.resize(m_noRkStages);
690 m_timeIntegrationCoefficientsC.resize(m_noRkStages);
692 m_timeIntegrationCoefficientsA = {{F0, 0.718801210867241, 0.778533117342157, 0.005328279665404, 0.855297993402928,
693 3.95641382457746, 1.57805753805874, 2.08370945525741, 0.748333418276161,
694 0.703286110656336, -0.001391709611768, 0.093207536963746, 0.951420047087595,
696 m_timeIntegrationCoefficientsB = {{0.036776245431967, 0.313629660755396, 0.153184869186903, 0.003009708681818,
697 0.332629379064611, 0.244025140535086, 0.371887923959228, 0.620412622158244,
698 0.152404317302874, 0.076089492741927, 0.007760421404098, 0.002464728475538,
699 0.078034834004939, 5.50597772702696}};
700 m_timeIntegrationCoefficientsC = {{F0, 0.036776245431967, 0.124968526272503, 0.244617770227770, 0.247614953107042,
701 0.296931112038247, 0.397814964580264, 0.527085458944033, 0.698126999417570,
702 0.819089083535213, 0.852705988709862, 0.860471181746282, 0.862706037696998,
710 m_timeIntegrationCoefficientsA.resize(m_noRkStages);
711 m_timeIntegrationCoefficientsB.resize(m_noRkStages);
712 m_timeIntegrationCoefficientsC.resize(m_noRkStages);
714 m_timeIntegrationCoefficientsA = {{F0, 0.6160178650170565, 0.4449487060774118, 1.0952033345276178,
715 1.2256030785959187, 0.2740182222332805, 0.0411952089052647, 0.1797084899153560,
716 1.1771530652064288, 0.4078831463120878, 0.8295636426191777, 4.7895970584252288,
717 0.6606671432964504}};
718 m_timeIntegrationCoefficientsB = {{0.0271990297818803, 0.1772488819905108, 0.0378528418949694, 0.6086431830142991,
719 0.2154313974316100, 0.2066152563885843, 0.0415864076069797, 0.0219891884310925,
720 0.9893081222650993, 0.0063199019859826, 0.3749640721105318, 1.6080235151003195,
721 0.0961209123818189}};
722 m_timeIntegrationCoefficientsC = {{F0, 0.0271990297818803, 0.0952594339119365, 0.1266450286591127,
723 0.1825883045699772, 0.3737511439063931, 0.5301279418422206, 0.5704177433952291,
724 0.5885784947099155, 0.6160769826246714, 0.6223252334314046, 0.6897593128753419,
725 0.9126827615920843}};
732 m_timeIntegrationCoefficientsA.resize(m_noRkStages);
733 m_timeIntegrationCoefficientsB.resize(m_noRkStages);
734 m_timeIntegrationCoefficientsC.resize(m_noRkStages);
736 m_timeIntegrationCoefficientsA = {{F0, 0.808316387498383, 1.503407858773331, 1.053064525050744, 1.463149119280508,
737 0.659288128108783, 1.667891931891068}};
738 m_timeIntegrationCoefficientsB = {{0.0119705267309784, 0.8886897793820711, 0.4578382089261419, 0.5790045253338471,
739 0.3160214638138484, 0.2483525368264122, 0.0677123095940884}};
740 m_timeIntegrationCoefficientsC = {{F0, 0.0119705267309784, 0.1823177940361990, 0.5082168062551849,
741 0.6532031220148590, 0.8534401385678250, 0.9980466084623790}};
748 m_timeIntegrationCoefficientsA.resize(m_noRkStages);
749 m_timeIntegrationCoefficientsB.resize(m_noRkStages);
750 m_timeIntegrationCoefficientsC.resize(m_noRkStages);
752 m_timeIntegrationCoefficientsA = {{F0, 0.5534431294501569, -0.0106598757020349, 0.5515812888932000,
753 1.8857903775587410, 5.7012957427932640, -2.1139039656647930,
754 0.5339578826675280}};
755 m_timeIntegrationCoefficientsB = {{0.0803793688273695, 0.5388497458569843, 0.0197497440903196, 0.0991184129733997,
756 0.7466920411064123, 1.6795842456188940, 0.2433728067008188,
757 0.1422730459001373}};
758 m_timeIntegrationCoefficientsC = {{F0, 0.0803793688273695, 0.3210064250338430, 0.3408501826604660,
759 0.3850364824285470, 0.5040052477534100, 0.6578977561168540,
760 0.9484087623348481}};
766 if(m_noRkStages == -1) {
767 TERMM(1,
"Invalid property value for time integration scheme");
770 MString dgPolynomialType =
"DG_POLY_LEGENDRE";
772 string2enum(Context::getSolverProperty<MString>(
"dgPolynomialType", m_solverId, AT_, &dgPolynomialType));
787 m_initPolyDeg = Context::getSolverProperty<MInt>(
"initPolyDeg", m_solverId, AT_, &m_initPolyDeg);
789 if(m_initPolyDeg < 0) {
790 mTerm(1, AT_,
"Negative polynomial order set in properties.");
808 m_minPolyDeg = numeric_limits<MInt>::max();
809 m_minPolyDeg = Context::getSolverProperty<MInt>(
"minPolyDeg", m_solverId, AT_, &m_minPolyDeg);
827 m_maxPolyDeg = Context::getSolverProperty<MInt>(
"maxPolyDeg", m_solverId, AT_, &m_maxPolyDeg);
832 m_initPolyDeg = max(1, m_initPolyDeg);
833 m_maxPolyDeg = max(1, m_maxPolyDeg);
834 m_minPolyDeg = max(1, m_minPolyDeg);
836 m_initPolyDeg = max(0, m_initPolyDeg);
837 m_maxPolyDeg = max(0, m_maxPolyDeg);
838 m_minPolyDeg = max(0, m_minPolyDeg);
841 m_minPolyDeg = min(m_minPolyDeg, m_initPolyDeg);
842 m_maxPolyDeg = max(m_maxPolyDeg, m_initPolyDeg);
846 m_initNoNodes1D = -1;
847 m_initNoNodes1D = Context::getSolverProperty<MInt>(
"initNoNodes", m_solverId, AT_, &m_initNoNodes1D);
849 m_minNoNodes1D = Context::getSolverProperty<MInt>(
"minNoNodes", m_solverId, AT_, &m_initNoNodes1D);
851 m_maxNoNodes1D = Context::getSolverProperty<MInt>(
"maxNoNodes", m_solverId, AT_, &m_initNoNodes1D);
854 m_initNoNodes1D = max(2, m_initNoNodes1D);
855 m_maxNoNodes1D = max(2, m_maxNoNodes1D);
856 m_minNoNodes1D = max(2, m_minNoNodes1D);
858 m_minNoNodes1D = min(m_minNoNodes1D, m_initNoNodes1D);
859 m_maxNoNodes1D = max(m_maxNoNodes1D, m_initNoNodes1D);
862 m_initNoNodes1D = m_initPolyDeg + 1;
863 m_minNoNodes1D = m_minPolyDeg + 1;
864 m_maxNoNodes1D = m_maxPolyDeg + 1;
880 m_calcErrorNorms = 1;
881 m_calcErrorNorms = Context::getSolverProperty<MBool>(
"calcErrorNorms", m_solverId, AT_, &m_calcErrorNorms);
883 if(m_calcErrorNorms) {
897 m_noAnalysisNodes = 2 * (m_maxPolyDeg + 1);
898 m_noAnalysisNodes = Context::getSolverProperty<MInt>(
"noAnalysisNodes", m_solverId, AT_, &m_noAnalysisNodes);
899 if(m_noAnalysisNodes < 2 * m_maxPolyDeg) {
901 "WARNING: Polynomial degree for error norms should be at least twice the maximum polynomial degree";
902 m_log << warning << std::endl;
903 cerr0 << warning << std::endl;
905 m_polyDegAnalysis = m_noAnalysisNodes - 1;
908 m_noAnalysisNodes = m_initNoNodes1D;
909 m_polyDegAnalysis = m_initPolyDeg;
925 m_noErrorDigits = Context::getSolverProperty<MInt>(
"noErrorDigits", m_solverId, AT_, &m_noErrorDigits);
936 m_useSponge = Context::getSolverProperty<MBool>(
"useSponge", m_solverId, AT_, &m_useSponge);
939 MString nameCoords =
"prefCoordinates_" + to_string(count);
940 MString namePolyDeg =
"prefPolyDeg_" + to_string(count);
941 MString nameOperator =
"prefOperator_" + to_string(count);
942 MString nameNoNodes1D =
"prefNoNodes_" + to_string(count);
947 const MInt polyDeg = Context::getSolverProperty<MInt>(namePolyDeg, m_solverId, AT_);
949 if(polyDeg < m_minPolyDeg || polyDeg > m_maxPolyDeg) {
950 TERMM(1,
"p-refinement: Polynomial degree of patch " + to_string(count)
951 +
" doesn't fit into range defined by minimum and maximum "
952 "polynomial degree");
958 TERMM(1,
"p-refinement: wrong number of coordinates for patch " + to_string(count));
961 array<MFloat, 2 * nDim> coords;
963 for(
MInt i = 0; i < 2 * nDim; i++) {
964 coords[i] = Context::getSolverProperty<MFloat>(nameCoords, m_solverId, AT_, i);
968 prefOperator = Context::getSolverProperty<MString>(nameOperator, m_solverId, AT_, &prefOperator);
970 MInt prefNoNodes1D = polyDeg + 1;
971 prefNoNodes1D = Context::getSolverProperty<MInt>(nameNoNodes1D, m_solverId, AT_, &prefNoNodes1D);
973 m_prefPatchesCoords.push_back(coords);
974 m_prefPatchesPolyDeg.push_back(polyDeg);
975 m_prefPatchesOperators.push_back(prefOperator);
976 m_prefPatchesNoNodes1D.push_back(prefNoNodes1D);
981 nameCoords =
"prefCoordinates_" + to_string(count);
982 namePolyDeg =
"prefPolyDeg_" + to_string(count);
983 nameOperator =
"prefOperator_" + to_string(count);
984 nameNoNodes1D =
"prefNoNodes_" + to_string(count);
998 m_pref = Context::getSolverProperty<MInt>(
"pref", m_solverId, AT_, &m_pref);
1007 MString adaptiveMethod =
"DG_ADAPTIVE_NONE";
1008 m_adaptivePref =
string2enum(Context::getSolverProperty<MString>(
"adaptivePref", m_solverId, AT_, &adaptiveMethod));
1010 if(hasAdaptivePref()) {
1012 TERMM(1,
"Multi-solver does not work with adaptive p-refinement yet");
1021 m_adaptiveInterval = Context::getSolverProperty<MInt>(
"adaptiveInterval", m_solverId, AT_, &m_adaptiveInterval);
1030 m_adaptiveThreshold =
1031 Context::getSolverProperty<MFloat>(
"adaptiveThreshold", m_solverId, AT_, &m_adaptiveThreshold);
1032 if(m_adaptiveInterval < 0 || m_adaptiveThreshold < 0) {
1033 stringstream errorMessage;
1034 errorMessage <<
"Error: Adaptive setting variables cannot be set to < 0 "
1035 "(adaptiveInterval/adaptiveThreshold)"
1037 TERMM(1, errorMessage.str());
1042 m_useCutOffBoundaries =
false;
1043 fill(m_cutOffBoundaryConditionIds.begin(), m_cutOffBoundaryConditionIds.end(), -1);
1060 m_useCutOffBoundaries =
true;
1061 m_useCutOffBoundaries =
1062 Context::getSolverProperty<MBool>(
"useCutOffBoundaries", m_solverId, AT_, &m_useCutOffBoundaries);
1065 if(m_useCutOffBoundaries) {
1066 for(
MInt i = 0; i < 2 * nDim; i++) {
1081 m_cutOffBoundaryConditionIds[i] = Context::getSolverProperty<MInt>(
"cutOffBoundaryConditionIds", m_solverId,
1082 AT_, &m_cutOffBoundaryConditionIds[i], i);
1096 m_useSourceRampUp =
false;
1097 m_useSourceRampUp = Context::getSolverProperty<MBool>(
"useSourceRampUp", solverId(), AT_, &m_useSourceRampUp);
1099 if(m_useSourceRampUp) {
1109 m_sourceRampUpTime = Context::getSolverProperty<MFloat>(
"sourceRampUpTime", solverId(), AT_);
1124template <MInt nDim,
class SysEqn>
1129 MInt noElements = 0;
1130 for(
MInt cellId = 0; cellId < grid().noInternalCells(); cellId++) {
1132 if(needElementForCell(cellId)) {
1138 if(m_maxNoSurfaces == -1) {
1140 m_maxNoSurfaces = calculateNeededNoSurfaces() + 100;
1141 m_log <<
"Max. no surfaces was not specified in the properties file (or "
1142 "set to -1). It was therefore calculated automatically for the "
1144 << m_maxNoSurfaces << endl;
1148 MInt maxNoElements = noElements;
1150 "MPI_IN_PLACE",
"maxNoElements");
1151 MInt minNoElements = noElements;
1153 "MPI_IN_PLACE",
"minNoElements");
1154 MInt maxNoSurfaces = m_maxNoSurfaces;
1156 "MPI_IN_PLACE",
"maxNoSurfaces");
1159 message <<
"Solver #" << solverId() <<
" - maximum number of DG elements among ranks: " << maxNoElements
1161 message <<
"Solver #" << solverId() <<
" - minimum number of DG elements among ranks: " << minNoElements
1163 message <<
"Solver #" << solverId() <<
" - maximum number of DG surfaces among ranks: " << maxNoSurfaces
1170 m_elements.maxPolyDeg(m_maxPolyDeg);
1171 m_elements.maxNoNodes1D(m_maxNoNodes1D);
1172 m_elements.noNodeVars(SysEqn::noNodeVars());
1173 m_elements.reset(noElements);
1176 m_surfaces.maxPolyDeg(m_maxPolyDeg);
1177 m_surfaces.maxNoNodes1D(m_maxNoNodes1D);
1178 m_surfaces.noNodeVars(SysEqn::noNodeVars());
1179 m_surfaces.reset(m_maxNoSurfaces);
1182 MInt noHElements = 0;
1183 for(
MInt cellId = 0; cellId < grid().noCells(); cellId++) {
1184 if(needHElementForCell(cellId)) {
1190 m_helements.reset(noHElements);
1193template <MInt nDim,
class SysEqn>
1195 vector<MFloat> L2Error(m_sysEqn.noVars());
1196 vector<MFloat> LInfError(m_sysEqn.noVars());
1197 vector<MFloat> L2ErrLocal(m_sysEqn.noVars());
1198 vector<MFloat> LInfErrLocal(m_sysEqn.noVars());
1199 calcErrorNorms(m_finalTime, L2Error, LInfError, L2ErrLocal, LInfErrLocal);
1201 cout <<
"WRITE EOC" << endl;
1203 if(m_calcErrorNorms && isMpiRoot()) {
1219template <MInt nDim,
class SysEqn>
1223 RECORD_TIMER_START(m_timers[Timers::Run]);
1224 RECORD_TIMER_START(m_timers[Timers::RunInit]);
1227 initSolverObjects();
1235 RECORD_TIMER_STOP(m_timers[Timers::RunInit]);
1238 MBool finalTimeStep =
false;
1239 while(!finalTimeStep) {
1240 finalTimeStep = step();
1246 RECORD_TIMER_STOP(m_timers[Timers::Run]);
1254template <MInt nDim,
class SysEqn>
1257 RECORD_TIMER_START(m_timers[Timers::InitData]);
1260 if(!m_isInitSolver) {
1261 TERMM(1,
"Solver was not initialized.");
1265 m_log <<
"Restart is enabled." << endl;
1271 if(SysEqn::noNodeVars() > 0 && !SysEqn::hasTimeDependentNodeVars()) {
1280 m_time = m_startTime;
1281 m_firstTimeStep =
true;
1288 if(SysEqn::noNodeVars() > 0) {
1289 updateNodeVariables();
1292 extendNodeVariables();
1296 m_isInitData =
true;
1298 RECORD_TIMER_STOP(m_timers[Timers::InitData]);
1306template <MInt nDim,
class SysEqn>
1309 RECORD_TIMER_START(m_timers[Timers::InitMainLoop]);
1313 TERMM(1,
"Solver data was not initialized.");
1318 if(m_writeInitialSolution) {
1323 if(SysEqn::noNodeVars() > 0 && !SysEqn::hasTimeDependentNodeVars()) {
1324 if(m_writeNodeVarsFile) {
1330 m_pointData.save(
false);
1331 m_surfaceData.save(
false);
1332 m_volumeData.save(
false);
1335 m_slice.save(
false);
1339 m_endAutoSaveTime = -1;
1340 m_endAutoSaveWritten =
false;
1341 const char* envJobEndTime = getenv(
"MAIA_JOB_END_TIME");
1345 const MString jobEndTime(envJobEndTime);
1346 MBool onlyDigits =
true;
1347 for(
auto&& character : jobEndTime) {
1348 if(!isdigit(character)) {
1349 m_log <<
"Warning: the environment variable MAIA_JOB_END_TIME "
1350 "included non-digit characters.\n"
1351 <<
"Saving a restart file before the compute job ends is NOT "
1360 m_endAutoSaveTime = stoi(jobEndTime);
1361 m_endAutoSaveTime -= m_noMinutesEndAutoSave * 60;
1362 m_log <<
"Activated automatic restart file writing at " << ctime(&m_endAutoSaveTime)
1363 <<
" (in Unix time: " << m_endAutoSaveTime <<
")." << endl;
1368 outputInitSummary();
1371 analyzeTimeStep(m_time, F0, F0, m_timeStep, F0);
1377 m_noAnalyzeTimeSteps = 0;
1380 m_isInitMainLoop =
true;
1382 RECORD_TIMER_STOP(m_timers[Timers::InitMainLoop]);
1387 RECORD_TIMER_START(m_timers[Timers::MainLoop]);
1388 RECORD_TIMER_START(m_timers[Timers::RungeKuttaStep]);
1389 RECORD_TIMER_START(m_timers[Timers::TimeDeriv]);
1390 prolongToSurfaces();
1391 applyForwardProjection();
1392 startMpiSurfaceExchange();
1393 RECORD_TIMER_STOP(m_timers[Timers::TimeDeriv]);
1394 RECORD_TIMER_STOP(m_timers[Timers::RungeKuttaStep]);
1395 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
1415template <MInt nDim,
class SysEqn>
1418 TERMM(1,
"deprecated");
1421 if(!m_isInitMainLoop) {
1422 TERMM(1,
"Main loop was not initialized.");
1425 RECORD_TIMER_START(m_timers[Timers::MainLoop]);
1427 startLoadTimer(AT_);
1430 if(m_rkStage == 0) {
1435 if(externalDt < 0.0) {
1437 if(m_firstTimeStep ==
true || (m_calcTimeStepInterval > 0 && m_timeStep % m_calcTimeStepInterval == 0)) {
1438 m_dt = calcTimeStep();
1439 m_firstTimeStep =
false;
1446 if(isAdaptationTimeStep(m_timeStep)) {
1448 cancelMpiRequests();
1449 adaptiveRefinement(m_timeStep);
1454 MBool finalTimeStep =
false;
1455 if(m_finalTime - m_time - m_dt < 1.0E-10) {
1456 finalTimeStep =
true;
1457 m_dt = m_finalTime - m_time;
1458 }
else if(m_timeStep == m_timeSteps) {
1459 finalTimeStep =
true;
1460 }
else if(m_timeStep == std::max(m_restartTimeStep, 0) + m_timeSteps) {
1463 finalTimeStep =
true;
1468 timeStepRk(m_time, m_dt);
1473 timeStepRk(m_time, m_dt, m_rkStage);
1475 m_rkStage = (m_rkStage + 1) % m_noRkStages;
1481 if(m_rkStage != 0) {
1482 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
1483 return finalTimeStep;
1491 m_noAnalyzeTimeSteps++;
1494 if((m_solutionInterval > 0 && m_timeStep % m_solutionInterval == 0) || (finalTimeStep && m_alwaysSaveFinalSolution)) {
1495 RECORD_TIMER_START(m_timers[Timers::MainLoopIO]);
1502 m_outputTime += tOutputEnd - tOutputStart;
1504 RECORD_TIMER_STOP(m_timers[Timers::MainLoopIO]);
1508 if((m_restartInterval > 0 && m_timeStep % m_restartInterval == 0) || (finalTimeStep && m_alwaysSaveFinalRestart)) {
1509 RECORD_TIMER_START(m_timers[Timers::MainLoopIO]);
1516 m_outputTime += tOutputEnd - tOutputStart;
1518 RECORD_TIMER_STOP(m_timers[Timers::MainLoopIO]);
1522 if(m_endAutoSaveTime != -1 && !m_endAutoSaveWritten && m_timeStep % m_endAutoSaveCheckInterval == 0) {
1524 MBool writeEndAutoSave =
false;
1526 && m_endAutoSaveTime
1527 <= chrono::duration_cast<chrono::seconds>(chrono::system_clock::now().time_since_epoch()).count()) {
1528 writeEndAutoSave =
true;
1530 MPI_Bcast(&writeEndAutoSave, 1, MPI_C_BOOL, 0, mpiComm(), AT_,
"writeEndAutoSave");
1532 if(writeEndAutoSave) {
1533 RECORD_TIMER_START(m_timers[Timers::MainLoopIO]);
1538 std::chrono::duration_cast<std::chrono::seconds>(std::chrono::system_clock::now().time_since_epoch()).count();
1539 m_log <<
"Finished end auto save at " << std::ctime(&now) <<
" (in Unix time: " << now <<
")." << endl;
1540 time_t jobEndTime = m_endAutoSaveTime + m_noMinutesEndAutoSave * 60;
1541 m_log <<
"Job will end at " << std::ctime(&jobEndTime) <<
" (in Unix time: " << jobEndTime <<
")." << endl;
1543 m_endAutoSaveWritten =
true;
1548 m_outputTime += tOutputEnd - tOutputStart;
1550 RECORD_TIMER_STOP(m_timers[Timers::MainLoopIO]);
1554 RECORD_TIMER_START(m_timers[Timers::MainLoopIO]);
1556 m_pointData.save(finalTimeStep);
1557 m_surfaceData.save(finalTimeStep);
1558 m_volumeData.save(finalTimeStep);
1560 m_slice.save(finalTimeStep);
1561 RECORD_TIMER_STOP(m_timers[Timers::MainLoopIO]);
1564 if((m_analysisInterval > 0 && m_timeStep % m_analysisInterval == 0) || (finalTimeStep)) {
1565 RECORD_TIMER_START(m_timers[Timers::Analysis]);
1567 const MFloat timeDiff =
wallTime() - m_analyzeTimeStart - m_outputTime;
1568 const MFloat runTimeRelative = timeDiff * noDomains() / m_noAnalyzeTimeSteps / m_statGlobalNoActiveDOFs;
1571 analyzeTimeStep(m_time, runTimeRelative, runTimeTotal, m_timeStep, m_dt);
1572 if(finalTimeStep && isMpiRoot()) {
1574 cout <<
"----------------------------------------"
1575 <<
"----------------------------------------" << endl;
1576 cout <<
" MAIA finished. Final time: " << m_time <<
" Time steps: " << m_timeStep << endl;
1577 cout <<
"----------------------------------------"
1578 <<
"----------------------------------------" << endl;
1584 m_noAnalyzeTimeSteps = 0;
1586 RECORD_TIMER_STOP(m_timers[Timers::Analysis]);
1598 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
1602 return finalTimeStep;
1607template <MInt nDim,
class SysEqn>
1610 RECORD_TIMER_START(m_timers[Timers::MainLoop]);
1613 if(!m_isInitMainLoop) {
1614 TERMM(1,
"Main loop was not initialized.");
1618 if(m_rkStage != 0) {
1620 "preTimeStep should only be called at the start of a new timestep: rkStage = " + std::to_string(m_rkStage));
1625 ASSERT(m_timeStep ==
globalTimeStep,
"Error: time step inconsistent.");
1628 if(m_firstTimeStep ==
true || (m_calcTimeStepInterval > 0 && m_timeStep % m_calcTimeStepInterval == 0)) {
1629 m_dt = calcTimeStep();
1630 m_firstTimeStep =
false;
1634 if(isAdaptationTimeStep(m_timeStep)) {
1636 cancelMpiRequests();
1638 adaptiveRefinement(m_timeStep);
1642 m_finalTimeStep =
false;
1643 if(m_finalTime - m_time - m_dt < 1.0E-10) {
1644 m_finalTimeStep =
true;
1645 m_dt = m_finalTime - m_time;
1646 }
else if(
g_multiSolverGrid && m_timeStep == std::max(m_restartTimeStep, 0) + m_timeSteps) {
1647 m_finalTimeStep =
true;
1650 m_finalTimeStep =
true;
1654 resetExternalSources();
1656 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
1661template <MInt nDim,
class SysEqn>
1664 RECORD_TIMER_START(m_timers[Timers::MainLoop]);
1667 timeStepRk(m_time, m_dt, m_rkStage);
1669 m_rkStage = (m_rkStage + 1) % m_noRkStages;
1671 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
1673 return (m_rkStage == 0);
1678template <MInt nDim,
class SysEqn>
1681 RECORD_TIMER_START(m_timers[Timers::MainLoop]);
1686 m_noAnalyzeTimeSteps++;
1689 if((m_analysisInterval > 0 && m_timeStep % m_analysisInterval == 0) || (m_finalTimeStep)) {
1690 RECORD_TIMER_START(m_timers[Timers::Analysis]);
1692 const MFloat timeDiff =
wallTime() - m_analyzeTimeStart - m_outputTime;
1693 const MFloat runTimeRelative = timeDiff * noDomains() / m_noAnalyzeTimeSteps / m_statGlobalNoActiveDOFs;
1696 analyzeTimeStep(m_time, runTimeRelative, runTimeTotal, m_timeStep, m_dt);
1697 if(m_finalTimeStep && isMpiRoot()) {
1699 cout <<
"----------------------------------------"
1700 <<
"----------------------------------------" << endl;
1701 cout <<
" MAIA finished. Final time: " << m_time <<
" Time steps: " << m_timeStep << endl;
1702 cout <<
"----------------------------------------"
1703 <<
"----------------------------------------" << endl;
1709 m_noAnalyzeTimeSteps = 0;
1711 RECORD_TIMER_STOP(m_timers[Timers::Analysis]);
1724 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
1729template <MInt nDim,
class SysEqn>
1733 if(!isActive())
return false;
1735 if((m_restartInterval > 0 && m_timeStep % m_restartInterval == 0) || (m_finalTimeStep && m_alwaysSaveFinalRestart)) {
1739 return writeRestart;
1744template <MInt nDim,
class SysEqn>
1746 const MBool NotUsed(writeBackup),
1747 const MString NotUsed(gridFileName),
1748 MInt* NotUsed(recalcIdTree)) {
1751 RECORD_TIMER_START(m_timers[Timers::MainLoop]);
1755 RECORD_TIMER_START(m_timers[Timers::MainLoopIO]);
1762 m_outputTime += tOutputEnd - tOutputStart;
1764 RECORD_TIMER_STOP(m_timers[Timers::MainLoopIO]);
1808 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
1813template <MInt nDim,
class SysEqn>
1817 RECORD_TIMER_START(m_timers[Timers::MainLoop]);
1818 RECORD_TIMER_START(m_timers[Timers::MainLoopIO]);
1821 if((m_solutionInterval > 0 && m_timeStep % m_solutionInterval == 0)
1822 || ((m_finalTimeStep || finalTimeStep) && m_alwaysSaveFinalSolution) || forceOutput) {
1828 m_outputTime += tOutputEnd - tOutputStart;
1832 m_pointData.save(finalTimeStep);
1833 m_surfaceData.save(finalTimeStep);
1834 m_volumeData.save(finalTimeStep);
1836 m_slice.save(finalTimeStep);
1838 RECORD_TIMER_STOP(m_timers[Timers::MainLoopIO]);
1839 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
1844template <MInt nDim,
class SysEqn>
1847 RECORD_TIMER_START(m_timers[Timers::Run]);
1852 TERMM_IF_COND(m_slice.enabled(),
"Fixme: DG slices with inactive ranks not supported.");
1856 RECORD_TIMER_START(m_timers[Timers::RunInit]);
1859 initSolverObjects();
1868 RECORD_TIMER_STOP(m_timers[Timers::RunInit]);
1873template <MInt nDim,
class SysEqn>
1878 if(!isActive())
return;
1880 RECORD_TIMER_START(m_timers[Timers::RunInit]);
1883 RECORD_TIMER_STOP(m_timers[Timers::RunInit]);
1894template <MInt nDim,
class SysEqn>
1897 RECORD_TIMER_START(m_timers[Timers::InitSolverObjects]);
1911 initPolynomialDegree();
1914 initInterpolation();
1917 initNodeCoordinates();
1923 checkCellProperties();
1929 if(hasMpiExchange()) {
1935 m_log <<
"Initializing sponge... ";
1936 sponge().init(m_maxPolyDeg, &grid(), &m_elements, &m_surfaces, &m_boundaryConditions, &m_sysEqn, mpiComm());
1937 m_log <<
"done" << endl;
1942 if(hasMpiExchange() && hasPref()) {
1943 exchangeMpiSurfacePolyDeg();
1947 initMortarProjections();
1952 m_surfaceData.init();
1954 m_volumeData.init();
1960 initSimulationStatistics();
1963 m_isInitSolver =
true;
1965 RECORD_TIMER_STOP(m_timers[Timers::InitSolverObjects]);
1973template <MInt nDim,
class SysEqn>
2009 if(!(Context::getSolverProperty<MBool>(
"loadGridMap", m_solverId, AT_))) {
2027 MString gridMapFileName = outputDir() +
"gridmap" + ParallelIo::fileExt();
2028 gridMapFileName = Context::getSolverProperty<MString>(
"gridMapFileName", m_solverId, AT_, &gridMapFileName);
2030 loadGridMap(gridMapFileName);
2047 check = Context::getSolverProperty<MBool>(
"checkGridMap", m_solverId, AT_, &check);
2050 const MString donorGridFileName = Context::getSolverProperty<MString>(
"donorGridFileName", m_solverId, AT_);
2052 checkGridMap(donorGridFileName);
2063template <MInt nDim,
class SysEqn>
2067 m_log <<
"Loading grip map from " << gridMapFileName <<
"..." << endl;
2070 const MInt noCells = grid().noInternalCells();
2073 using namespace parallel_io;
2074 ParallelIo gridMapFile(gridMapFileName, PIO_READ, mpiComm());
2075 gridMapFile.setOffset(noCells, grid().domainOffset(domainId()));
2076 gridMapFile.readArray(&gridMap[0],
"cellIds");
2077 gridMapFile.readArray(&gridMapNoOffspring[0],
"noOffspring");
2083 m_gridMapOffsets.clear();
2086 for(
MInt cellId = 0; cellId < noCells; cellId++) {
2088 if(gridMap[cellId] == -1) {
2093 const MInt firstTargetCellId = cellId;
2094 const MInt firstDonorGlobalId = gridMap[cellId];
2098 noOffspring[0] = gridMapNoOffspring[cellId];
2103 MInt noDonorCells = 1 + gridMapNoOffspring[cellId];
2104 while(cellId + 1 < noCells) {
2105 if(gridMap[cellId + 1] == firstDonorGlobalId + noDonorCells) {
2109 noDonorCells += 1 + gridMapNoOffspring[cellId];
2110 noOffspring[cellId - firstTargetCellId] = gridMapNoOffspring[cellId];
2111 }
else if(gridMap[cellId + 1] == firstDonorGlobalId + noDonorCells - 1) {
2116 noOffspring[cellId - firstTargetCellId] = -1;
2124 const MInt noTargetCells = cellId - firstTargetCellId + 1;
2127 m_gridMapOffsets.push_back({firstTargetCellId, noTargetCells, firstDonorGlobalId, noDonorCells,
2128 std::vector<MInt>(noOffspring.
begin(), noOffspring.
end())});
2131 m_log <<
"Found " << m_gridMapOffsets.size() <<
" contiguous mapped region(s):" << endl;
2132 for(
size_t i = 0; i < m_gridMapOffsets.size(); i++) {
2133 const auto& o = m_gridMapOffsets[i];
2134 m_log <<
"- offset: " << i <<
"; firstTargetCellId: " << o.m_firstTargetCellId
2135 <<
"; globalId: " << grid().tree().globalId(o.m_firstTargetCellId) <<
"; noDonorCells: " << o.m_noDonorCells
2136 <<
"; noTargetCells: " << o.m_noTargetCells <<
"; firstDonorGlobalId: " << o.m_firstDonorGlobalId << endl;
2140 m_maxNoGridMapOffsets = m_gridMapOffsets.size();
2142 "MPI_IN_PLACE",
"m_maxNoGridMapOffsets");
2143 m_log <<
"Maximum number of grid map offsets on all domains: " << m_maxNoGridMapOffsets << endl;
2156template <MInt nDim,
class SysEqn>
2163 if(domainId() == 0) {
2164 cerr <<
"Warning: DgCartesianSolver::checkGridMap deactivated! (fix coordinate comparison)" << endl;
2171 m_log <<
"Checking grip map for donor grid " << donorGridFileName <<
"...";
2176 MInt noDonorCellsLocal = 0;
2177 MInt gridMapMaxNoDonorCells = 0;
2178 for(
auto&& offset : m_gridMapOffsets) {
2179 noDonorCellsLocal += offset.m_noDonorCells;
2180 gridMapMaxNoDonorCells = max(gridMapMaxNoDonorCells, offset.m_noDonorCells);
2188 ParallelIo file(donorGridFileName, PIO_READ, mpiComm());
2189 MInt noDonorCellsFile;
2191 file.getAttribute(&noDonorCellsFile,
"noCells");
2192 MInt noDonorCellsGlobal;
2194 "noDonorCellsLocal",
"noDonorCellsGlobal");
2195 if(noDonorCellsFile != noDonorCellsGlobal) {
2196 TERMM(1,
"Number of mapped donor cells does not match number of cells in "
2197 "donor grid file: mapped: "
2198 + to_string(noDonorCellsGlobal) +
"; grid: " + to_string(noDonorCellsFile));
2206 const MFloat eps = 1.0 /
FPOW2(30) * grid().lengthLevel0();
2211 const array<MString, 3> dirs = {{
"x",
"y",
"z"}};
2212 for(
MInt offsetId = 0; offsetId < m_maxNoGridMapOffsets; offsetId++) {
2215 const MBool valid = (
static_cast<size_t>(offsetId) < m_gridMapOffsets.size());
2218 const MInt noDonorCells = valid ? m_gridMapOffsets[offsetId].m_noDonorCells : 0;
2219 const MInt firstDonorGlobalId = valid ? m_gridMapOffsets[offsetId].m_firstDonorGlobalId : 0;
2220 file.setOffset(noDonorCells, firstDonorGlobalId);
2223 for(
MInt dir = 0; dir < nDim; dir++) {
2225 file.readArray(&coordinates[0],
"coordinates_" + to_string(dir));
2232 const MInt noTargetCells = m_gridMapOffsets[offsetId].m_noTargetCells;
2233 MInt donorCellId = 0;
2236 MFloat donorLength = -numeric_limits<MFloat>::infinity();
2239 for(
MInt i = 0; i < noTargetCells; i++) {
2240 const MInt cellId = m_gridMapOffsets[offsetId].m_firstTargetCellId + i;
2241 const MFloat targetCoordinate = grid().tree().coordinate(cellId, dir);
2244 if(m_gridMapOffsets[offsetId].m_noOffspring[i] == -1) {
2245 const MFloat lastDonorCoordinate = coordinates[donorCellId - 1];
2249 if(targetCoordinate < lastDonorCoordinate - 0.5 * donorLength
2250 || targetCoordinate > lastDonorCoordinate + 0.5 * donorLength) {
2251 TERMM(1,
"Target cell of one-to-many mapping not contained in "
2259 const MFloat donorCoordinate = coordinates[donorCellId];
2260 if(!
approx(targetCoordinate, donorCoordinate, eps)) {
2261 const MInt globalId = grid().tree().globalId(cellId);
2262 TERMM(1, dirs[dir] +
"-coordinates for target cell " + to_string(cellId) +
" (globalId: "
2263 + to_string(globalId) +
") and donor cell " + to_string(firstDonorGlobalId + donorCellId)
2264 +
" do not match: target: " + to_string(targetCoordinate)
2265 +
"; donor: " + to_string(donorCoordinate));
2269 const MFloat targetLength = grid().cellLengthAtCell(cellId);
2271 donorLength = targetLength;
2274 for(
MInt offspringId = 0; offspringId < m_gridMapOffsets[offsetId].m_noOffspring[i]; offspringId++) {
2275 const MFloat offspringCoordinate = coordinates[donorCellId + offspringId];
2276 if(offspringCoordinate < targetCoordinate - 0.5 * targetLength
2277 || offspringCoordinate > targetCoordinate + 0.5 * targetLength) {
2278 TERMM(1,
"Donor offspring not contained in target cell.");
2282 donorCellId += m_gridMapOffsets[offsetId].m_noOffspring[i];
2288 m_log <<
"OK (duration: " << (
wallTime() - startTime) <<
" s)" << endl;
2297template <MInt nDim,
class SysEqn>
2301 m_log <<
"Initializing elements... ";
2302 for(
MInt cellId = 0; cellId < grid().noInternalCells(); cellId++) {
2304 if(needElementForCell(cellId)) {
2305 createElement(cellId);
2308 m_log <<
"done" << endl;
2316template <MInt nDim,
class SysEqn>
2320 m_log <<
"Initializing h-refined elements... ";
2322 for(
MInt cellId = 0; cellId < grid().noInternalCells(); cellId++) {
2324 if(needHElementForCell(cellId)) {
2325 createHElement(cellId);
2329 m_log <<
"done" << endl;
2341template <MInt nDim,
class SysEqn>
2345 MBool needElement =
true;
2349 if(grid().tree().hasChildren(cellId)
2351 || !pointIsInside(&grid().tree().coordinate(cellId, 0))) {
2352 needElement =
false;
2366template <MInt nDim,
class SysEqn>
2371 if(grid().tree().hasChildren(cellId)) {
2377 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
2378 if(grid().tree().hasNeighbor(cellId, dir)) {
2379 const MInt neighborId = grid().tree().neighbor(cellId, dir);
2380 if(grid().tree().hasChildren(neighborId)) {
2397template <MInt nDim,
class SysEqn>
2402 MInt noSurfaces = 0;
2403 for(
MInt cellId = 0; cellId < grid().noInternalCells(); cellId++) {
2405 if(grid().tree().hasChildren(cellId)) {
2410 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
2411 if(!grid().tree().hasNeighbor(cellId, dir)) {
2417 const MInt neighborId = grid().tree().neighbor(cellId, dir);
2418 const MBool neighborIsHalo = (neighborId >= grid().noInternalCells());
2420 if(grid().tree().hasChildren(neighborId)) {
2422 if(neighborIsHalo) {
2427 noSurfaces +=
IPOW2(nDim - 1);
2431 if(dir % 2 == 1 || neighborIsHalo) {
2459template <MInt nDim,
class SysEqn>
2464 const MFloat*
const box = geometry().boundingBox();
2467 MBool isInside =
true;
2468 for(
MInt dir = 0; dir < nDim; dir++) {
2470 array<MFloat, 2 * nDim> target;
2471 for(
MInt i = 0; i < nDim; i++) {
2472 target[i] = coordinates[i];
2473 target[i + nDim] = coordinates[i];
2476 target[dir] = box[dir] - (target[dir] - box[dir]);
2479 std::vector<MInt> nodeList;
2480 geometry().getLineIntersectionElements(&target[0], nodeList);
2483 if(nodeList.size() % 2 == 0) {
2499template <MInt nDim,
class SysEqn>
2504 const MInt elementId = m_elements.size();
2505 m_elements.append();
2508 m_elements.cellId(elementId) = cellId;
2511 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
2512 m_elements.surfaceIds(elementId, dir) = -1;
2522template <MInt nDim,
class SysEqn>
2527 const MInt hElementId = m_helements.size();
2528 m_helements.append();
2531 m_helements.elementId(hElementId) = m_elements.getElementByCellId(cellId);
2534 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
2535 for(
MInt pos = 0; pos < 2 * (nDim - 1); pos++) {
2536 m_helements.hrefSurfaceIds(hElementId, dir, pos) = -1;
2548template <MInt nDim,
class SysEqn>
2552 m_log <<
"Initializing polynomial degree and number of nodes in the elements... ";
2555 fill_n(&m_elements.polyDeg(0), m_elements.size(), m_initPolyDeg);
2556 fill_n(&m_elements.noNodes1D(0), m_elements.size(), m_initNoNodes1D);
2566 m_internalDataSize = m_elements.size() * m_elements.maxNoNodesXD() * SysEqn::noVars();
2571 m_noTotalNodesXD = 0;
2572 for(
MInt i = 0; i < m_elements.size(); i++) {
2573 m_noTotalNodesXD += m_elements.noNodesXD(i);
2576 m_log <<
"done" << endl;
2584template <MInt nDim,
class SysEqn>
2588 m_log <<
"Static p-refinement is enabled." << endl;
2591 fill(changed.
begin(), changed.
end(), 0);
2592 const MInt noElements = m_elements.size();
2594 for(
MInt elementId = 0; elementId < noElements; elementId++) {
2595 const MInt cellId = m_elements.cellId(elementId);
2603 MInt elemPolyDeg = m_initPolyDeg;
2604 MInt elemNoNodes1D = m_initNoNodes1D;
2608 for(std::vector<MFloat>::size_type patchId = 0; patchId < m_prefPatchesPolyDeg.size(); patchId++) {
2609 MBool inPatch =
true;
2611 for(
MInt i = 0; i < nDim; i++) {
2612 const MFloat cellCoord = grid().tree().coordinate(cellId, i);
2614 if(cellCoord < m_prefPatchesCoords[patchId][i]) {
2617 if(cellCoord > m_prefPatchesCoords[patchId][i + nDim]) {
2624 elemPolyDeg = m_prefPatchesPolyDeg[patchId];
2625 elemNoNodes1D = m_prefPatchesNoNodes1D[patchId];
2627 changed[elementId]++;
2630 m_elements.polyDeg(elementId) = elemPolyDeg;
2631 m_elements.noNodes1D(elementId) = elemNoNodes1D;
2634 m_log <<
"p-refinement"
2635 <<
" changed the polynomial degree "
2636 <<
"on " << count_if(changed.
begin(), changed.
end(), [](
MInt i) { return i > 0; }) <<
" elements in total for "
2637 << accumulate(changed.
begin(), changed.
end(), 0) <<
" times." << endl;
2648template <MInt nDim,
class SysEqn>
2652 m_log <<
"Initializing interpolation data... ";
2655 m_interpolation.clear();
2656 m_interpolation.resize(m_maxPolyDeg + 1);
2657 for(
MInt polyDeg = m_minPolyDeg; polyDeg <= m_maxPolyDeg; polyDeg++) {
2658 m_interpolation[polyDeg].clear();
2659 m_interpolation[polyDeg].resize(m_maxNoNodes1D + 1);
2669 for(
MInt i = m_minPolyDeg; i <= m_maxPolyDeg; i++) {
2670 MInt prefIndex = -1;
2671 for(std::vector<MFloat>::size_type j = 0; j < m_prefPatchesPolyDeg.size(); j++) {
2672 if(i == (
MInt)m_prefPatchesPolyDeg[j]) {
2678 MString refOperator = m_sbpOperator;
2679 MInt refNoNodes1D = m_sbpMode ? m_initNoNodes1D : (i + 1);
2680 if(prefIndex != -1) {
2681 refOperator = m_prefPatchesOperators[prefIndex];
2682 refNoNodes1D = m_prefPatchesNoNodes1D[prefIndex];
2684 m_interpolation[i][refNoNodes1D].init(i, polyType, refNoNodes1D, intMethod, m_sbpMode, refOperator);
2689 m_globalVolume = F0;
2691 for(
MInt elementId = 0; elementId < m_elements.size(); elementId++) {
2692 const MInt cellId = m_elements.cellId(elementId);
2693 m_localVolume += pow(grid().cellLengthAtCell(cellId), nDim);
2696 RECORD_TIMER_START(m_timers[Timers::Accumulated]);
2697 RECORD_TIMER_START(m_timers[Timers::MPI]);
2698 RECORD_TIMER_START(m_timers[Timers::MPIComm]);
2699 MPI_Allreduce(&m_localVolume, &m_globalVolume, 1, MPI_DOUBLE, MPI_SUM, mpiComm(), AT_,
"m_localVolume",
2701 RECORD_TIMER_STOP(m_timers[Timers::MPIComm]);
2702 RECORD_TIMER_STOP(m_timers[Timers::MPI]);
2703 RECORD_TIMER_STOP(m_timers[Timers::Accumulated]);
2706 if(m_calcErrorNorms) {
2708 m_interpAnalysis.init(m_polyDegAnalysis, polyType, m_noAnalysisNodes, intMethod, m_sbpMode, m_sbpOperator);
2711 const MInt noNodesAnalysis1D = m_noAnalysisNodes;
2712 const MInt noNodesAnalysis1D3 = (nDim == 3) ? noNodesAnalysis1D : 1;
2713 m_wVolumeAnalysis.resize(noNodesAnalysis1D, noNodesAnalysis1D, noNodesAnalysis1D3);
2715 for(
MInt i = 0; i < noNodesAnalysis1D; i++) {
2716 for(
MInt j = 0; j < noNodesAnalysis1D; j++) {
2717 for(
MInt k = 0; k < noNodesAnalysis1D3; k++) {
2718 m_wVolumeAnalysis(i, j, k) = wInt[i] * wInt[j] * (nDim == 3 ? wInt[k] : F1);
2725 m_vdmAnalysis.
clear();
2726 m_vdmAnalysis.resize(m_maxPolyDeg + 1);
2727 for(
MInt polyDeg = m_minPolyDeg; polyDeg <= m_maxPolyDeg; polyDeg++) {
2728 m_vdmAnalysis[polyDeg].clear();
2729 m_vdmAnalysis[polyDeg].resize(m_maxNoNodes1D + 1);
2734 for(
MInt polyDeg = m_minPolyDeg; polyDeg <= m_maxPolyDeg; polyDeg++) {
2736 for(
MInt noNodes1D = m_minNoNodes1D; noNodes1D <= m_maxNoNodes1D; noNodes1D++) {
2737 m_vdmAnalysis[polyDeg][noNodes1D].resize(m_noAnalysisNodes, noNodes1D);
2739 ASSERT(noNodes1D == m_noAnalysisNodes,
2740 "For now SBP Error Analysis only supports direct evaluation at regular nodes "
2741 "without interpolation...(noNodes1D = "
2742 + to_string(noNodes1D) +
", noAnalysisNodes = " + to_string(m_noAnalysisNodes) +
")");
2747 &m_interpAnalysis.m_nodes[0],
2748 &m_vdmAnalysis[polyDeg][noNodes1D][0]);
2751 const MInt noNodes1D = polyDeg + 1;
2752 m_vdmAnalysis[polyDeg][noNodes1D].resize(m_noAnalysisNodes, noNodes1D);
2757 &m_interpAnalysis.m_nodes[0],
2759 &m_vdmAnalysis[polyDeg][noNodes1D][0]);
2764 m_log <<
"done" << endl;
2774template <MInt nDim,
class SysEqn>
2778 m_log <<
"Initializing node coordinates... ";
2780 const MInt noElements = m_elements.size();
2783 for(
MInt elementId = 0; elementId < noElements; elementId++) {
2784 calcElementNodeCoordinates(elementId);
2787 m_log <<
"done" << endl;
2798template <MInt nDim,
class SysEqn>
2801 const MInt cellId = m_elements.cellId(elementId);
2802 const MInt lvl = grid().tree().level(cellId);
2803 const MInt polyDeg = m_elements.polyDeg(elementId);
2804 const MInt noNodes1D = m_elements.noNodes1D(elementId);
2805 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
2806 MFloatTensor nodeCoordinates(&m_elements.nodeCoordinates(elementId), noNodes1D, noNodes1D, noNodes1D3, nDim);
2809 for(
MInt i = 0; i < noNodes1D; i++) {
2810 for(
MInt j = 0; j < noNodes1D; j++) {
2811 for(
MInt k = 0; k < noNodes1D3; k++) {
2814 nodeCoordinates(i, j, k, 0) =
2815 grid().tree().coordinate(cellId, 0)
2816 + F1B2 * grid().cellLengthAtLevel(lvl) * m_interpolation[polyDeg][noNodes1D].m_nodes[i];
2817 nodeCoordinates(i, j, k, 1) =
2818 grid().tree().coordinate(cellId, 1)
2819 + F1B2 * grid().cellLengthAtLevel(lvl) * m_interpolation[polyDeg][noNodes1D].m_nodes[j];
2820 IF_CONSTEXPR(nDim == 3) {
2821 nodeCoordinates(i, j, k, 2) =
2822 grid().tree().coordinate(cellId, 2)
2823 + F1B2 * grid().cellLengthAtLevel(lvl) * m_interpolation[polyDeg][noNodes1D].m_nodes[k];
2844template <MInt nDim,
class SysEqn>
2848 m_log <<
"Initializing inverse Jacobian determinant... ";
2850 const MInt noElements = m_elements.size();
2851 MFloat* invJacobians = &m_elements.invJacobian(0);
2853 for(
MInt i = 0; i < noElements; i++) {
2854 const MInt cellId = m_elements.cellId(i);
2855 invJacobians[i] = F2 / grid().cellLengthAtCell(cellId);
2858 m_log <<
"done" << endl;
2869template <MInt nDim,
class SysEqn>
2873 m_log <<
"Initializing cell properties... ";
2875 const MInt noCells = grid().noCells();
2879 for(
MInt c = grid().noInternalCells(); c < noCells; c++) {
2880 if(!grid().tree().hasProperty(c, Cell::IsHalo)) {
2881 TERMM(1,
"Cell " + to_string(c) +
" should be marked as a halo cell");
2885 m_log <<
"done" << endl;
2892template <MInt nDim,
class SysEqn>
2896 const MInt noElements = m_elements.size();
2897 const MInt noDirs = 2 * nDim;
2898 const MInt*
const cellIds = &m_elements.cellId(0);
2903 m_log <<
"Creating boundary surfaces... ";
2905 m_noBoundarySurfaces = 0;
2910 this->identifyBoundaryCells(&isInterface[0]);
2913 if(m_geometryIntersection) {
2914 delete m_geometryIntersection;
2922 fill(bcIds.
begin(), bcIds.
end(), -1);
2923 for(
MInt elementId = 0; elementId < noElements; elementId++) {
2925 const MInt cellId = cellIds[elementId];
2926 if(!isInterface[cellId]) {
2931 array<MInt, 2 * nDim> ids = getBoundaryConditionIds(cellId);
2934 for(
MInt dir = 0; dir < noDirs; dir++) {
2935 const MInt nId = grid().tree().neighbor(cellId, dir);
2936 if(ids[dir] > -1 && nId > -1 && needElementForCell(nId)) {
2941 copy(ids.begin(), ids.end(), &bcIds(elementId, 0));
2949 for(
MInt elementId = 0; elementId < noElements; elementId++) {
2950 const MInt cellId = cellIds[elementId];
2951 for(
MInt dir = 0; dir < noDirs; dir++) {
2953 if(bcIds(elementId, dir) > -1) {
2958 if(!grid().tree().hasNeighbor(cellId, dir)) {
2963 const MInt neighborCellId = grid().tree().neighbor(cellId, dir);
2964 if(!isInterface[neighborCellId]) {
2969 if(needElementForCell(neighborCellId)) {
2974 auto ids = getBoundaryConditionIds(neighborCellId);
2978 const MInt oppositeDir = 2 * (dir / 2) + 1 - (dir % 2);
2979 if(ids[oppositeDir] > -1) {
2980 bcIds(elementId, dir) = ids[oppositeDir];
2987 if(m_useCutOffBoundaries) {
2989 fill(m_noCutOffBoundarySurfaces.begin(), m_noCutOffBoundarySurfaces.end(), 0);
2992 for(
MInt elementId = 0; elementId < noElements; elementId++) {
2993 const MInt cellId = cellIds[elementId];
2994 for(
MInt dir = 0; dir < noDirs; dir++) {
2996 if(m_cutOffBoundaryConditionIds[dir] < 0) {
3001 if(bcIds(elementId, dir) > -1) {
3006 if(grid().tree().hasNeighbor(cellId, dir)) {
3012 if(grid().tree().parent(cellId) > -1 && grid().tree().hasNeighbor(grid().tree().parent(cellId), dir)) {
3017 bcIds(elementId, dir) = m_cutOffBoundaryConditionIds[dir];
3020 m_noCutOffBoundarySurfaces[dir]++;
3026 set<MInt> localUniqueBcIds(bcIds.
begin(), bcIds.
end());
3029 localUniqueBcIds.erase(-1);
3037 MInt noBcIds = localUniqueBcIds.size();
3044 MPI_Allgather(&noBcIds, 1, MPI_INT, &countsBcIds[0], 1, MPI_INT, mpiComm(), AT_,
"noBcIds",
"countsBcIds[0]");
3047 offsetsBcIds[0] = 0;
3048 for(
MInt dId = 1; dId < noDomains(); dId++) {
3049 offsetsBcIds[dId] = offsetsBcIds[dId - 1] + countsBcIds[dId - 1];
3055 std::copy(localUniqueBcIds.begin(), localUniqueBcIds.end(), &localBcIds[0]);
3058 const MInt totalNoBcIds = offsetsBcIds[noDomains() - 1] + countsBcIds[noDomains() - 1];
3062 MPI_Allgatherv(&localBcIds[0], noBcIds, MPI_INT, &globalBcIds[0], &countsBcIds[0], &offsetsBcIds[0], MPI_INT,
3063 mpiComm(), AT_,
"localBcIds[0]",
"globalBcIds[0]");
3066 set<MInt> uniqueBcIds(globalBcIds.
begin(), globalBcIds.
end());
3070 std::vector<MInt> boundaryConditionIds;
3071 std::vector<std::pair<MInt, MInt>> bcIntervals;
3072 for(
const auto& boundaryConditionId : uniqueBcIds) {
3073 const MInt begin = m_noBoundarySurfaces;
3075 for(
MInt elementId = 0; elementId < noElements; elementId++) {
3076 for(
MInt dir = 0; dir < noDirs; dir++) {
3078 if(boundaryConditionId != bcIds(elementId, dir)) {
3081 if(hasSurface(elementId, dir)) {
3088 if(grid().periodicCartesianDir(dir / 2) != 0) {
3091 if(boundaryConditionId != 0) {
3093 "If you use periodic BCs, you must set the BC id of corresponding boundaries "
3099 initBoundarySurface(elementId, dir);
3102 const MInt end = m_noBoundarySurfaces;
3107 boundaryConditionIds.push_back(boundaryConditionId);
3108 bcIntervals.emplace_back(begin, end);
3111 m_log <<
"done" << endl;
3116 m_log <<
"creating inner surfaces... ";
3118 m_innerSurfacesOffset = m_surfaces.size();
3119 m_noInnerSurfaces = 0;
3121 for(
MInt elementId = 0; elementId < noElements; elementId++) {
3122 for(
MInt dir = 0; dir < noDirs; dir++) {
3123 if(hasSurface(elementId, dir)) {
3126 initInnerSurface(elementId, dir);
3130 m_log <<
"done" << endl;
3135 m_log <<
"creating MPI surfaces... ";
3137 m_mpiSurfacesOffset = m_surfaces.size();
3138 m_noMpiSurfaces = 0;
3140 if(hasMpiExchange()) {
3141 for(
MInt elementId = 0; elementId < noElements; elementId++) {
3142 for(
MInt dir = 0; dir < noDirs; dir++) {
3143 const MInt cellId = m_elements.cellId(elementId);
3144 const MInt nghbrId = grid().tree().neighbor(cellId, dir);
3145 MBool partLvlAncestorNghbr =
false;
3150 if(nghbrId > -1 && grid().tree().hasProperty(nghbrId, Cell::IsPartLvlAncestor)
3151 && grid().tree().hasChildren(nghbrId)) {
3154 for(
MInt child = 0; child <
IPOW2(nDim); child++) {
3155 const MInt childId = grid().tree().child(nghbrId, child);
3156 if(childId > -1 && grid().tree().hasProperty(childId, Cell::IsHalo)) {
3157 partLvlAncestorNghbr =
true;
3162 if(hasSurface(elementId, dir) && !partLvlAncestorNghbr) {
3165 initMpiSurface(elementId, dir);
3170 m_log <<
"done" << endl;
3173 m_boundaryConditions.clear();
3176 m_log <<
"creating and initializing boundary conditions... ";
3177 for(
MUint bcId = 0; bcId < boundaryConditionIds.size(); bcId++) {
3178 m_boundaryConditions.emplace_back(make_bc(boundaryConditionIds[bcId]));
3180 const MInt begin = bcIntervals[bcId].first;
3181 const MInt end = bcIntervals[bcId].second;
3182 m_boundaryConditions[bcId]->init(begin, end);
3184 m_log <<
"done" << endl;
3187 m_log <<
"Sanity check: all elements must have surfaces in each "
3190 for(
MInt elementId = 0; elementId < noElements; elementId++) {
3191 for(
MInt dir = 0; dir < noDirs; dir++) {
3192 if(!hasSurface(elementId, dir)) {
3194 m_log <<
"\nmissing surface for element " << elementId <<
" in direction " << dir <<
" (coordinates: ";
3195 const MInt cellId = m_elements.cellId(elementId);
3196 for(
MInt i = 0; i < nDim; i++) {
3197 m_log << grid().tree().coordinate(cellId, i) <<
" ";
3199 m_log <<
")" << std::endl;
3210 TERMM(1,
"There are " + to_string(missing)
3211 +
" surfaces missing (check m_log for more details). Did "
3212 "you maybe forget to specify cut-off boundary conditions?");
3214 m_log <<
"done" << endl;
3218 if(grid().periodicCartesianDir(0) && grid().periodicCartesianDir(1) && (nDim == 2 || grid().periodicCartesianDir(2))
3219 && m_noBoundarySurfaces > 0) {
3220 TERMM(1,
"There are boundary surfaces although every direction has a periodic boundaries. "
3221 "Please check what is wrong here.\n");
3235template <MInt nDim,
class SysEqn>
3238 return m_elements.surfaceIds(elementId, dir) > -1;
3250template <MInt nDim,
class SysEqn>
3254 const MInt surfaceId = createSurface(elementId, dir);
3255 m_noBoundarySurfaces++;
3278template <MInt nDim,
class SysEqn>
3282 MInt cellId = m_elements.cellId(elementId);
3285 if(!grid().tree().hasNeighbor(cellId, dir)) {
3286 const MInt parentId = grid().tree().parent(cellId);
3287 if(parentId > -1 && grid().tree().hasNeighbor(parentId, dir)) {
3298 const MInt nghbrId = grid().tree().neighbor(cellId, dir);
3301 if(grid().tree().hasChildren(nghbrId)) {
3306 if(grid().tree().hasProperty(nghbrId, Cell::IsHalo)) {
3311 const MInt nghbrElementId = m_elements.getElementByCellId(nghbrId);
3312 if(nghbrElementId == -1) {
3317 const MInt surfaceId = createSurface(elementId, dir);
3318 m_noInnerSurfaces++;
3337template <MInt nDim,
class SysEqn>
3341 MInt cellId = m_elements.cellId(elementId);
3344 if(!grid().tree().hasNeighbor(cellId, dir)) {
3345 const MInt parentId = grid().tree().parent(cellId);
3346 if(parentId > -1 && grid().tree().hasNeighbor(parentId, dir)) {
3357 const MInt nghbrId = grid().tree().neighbor(cellId, dir);
3362 if(!grid().tree().hasProperty(nghbrId, Cell::IsHalo) && !grid().tree().hasChildren(nghbrId)) {
3366 MInt surfaceId = -1;
3369 if(grid().tree().hasChildren(nghbrId)) {
3371 surfaceId = createHMPISurfaces(elementId, dir);
3374 surfaceId = createSurface(elementId, dir);
3394template <MInt nDim,
class SysEqn>
3400 array<MFloat, 2 * nDim> targetRegion;
3401 for(
MInt dir = 0; dir < nDim; dir++) {
3402 targetRegion[dir] = grid().tree().coordinate(cellId, dir) - 0.5 * grid().cellLengthAtCell(cellId);
3403 targetRegion[dir + nDim] = grid().tree().coordinate(cellId, dir) + 0.5 * grid().cellLengthAtCell(cellId);
3407 const MFloat cellHalfLength = 0.5 * grid().cellLengthAtCell(cellId);
3408 array<MFloat, nDim> coordinates;
3409 for(
MInt i = 0; i < nDim; i++) {
3410 coordinates[i] = grid().tree().coordinate(cellId, i);
3414 const MBool hasElement = needElementForCell(cellId);
3417 std::vector<MInt> nodeList;
3418 geometry().getIntersectionElements(&targetRegion[0], nodeList, cellHalfLength, &coordinates[0]);
3420 array<MInt, 2 * nDim> bcIds;
3424 static constexpr MFloat standardBasis[MAX_SPACE_DIMENSIONS][MAX_SPACE_DIMENSIONS] = {
3425 {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
3426 for(
MInt n = 0; n < (signed)nodeList.size(); n++) {
3430 array<MFloat, nDim * nDim> geometryPoints;
3434 array<MFloat, nDim> normal;
3435 elem.
calcNormal(&geometryPoints[0], &normal[0]);
3438 MInt orientation = -1;
3439 for(
MInt d = 0; d < nDim; d++) {
3440 const MFloat dot = inner_product(&normal[0], &normal[nDim], standardBasis[d], 0.0);
3441 if(
approx(std::fabs(dot), 1.0, MFloatEps)) {
3448 if(orientation != -1) {
3453 array<MFloat, nDim> centroid;
3455 const MFloat pos = grid().tree().coordinate(cellId, orientation) + 0.5 * grid().cellLengthAtCell(cellId);
3456 const MFloat neg = grid().tree().coordinate(cellId, orientation) - 0.5 * grid().cellLengthAtCell(cellId);
3457 const MInt direction = (std::fabs(neg - centroid[orientation]) < std::fabs(pos - centroid[orientation]))
3459 : 2 * orientation + 1;
3463 if(bcIds[direction] > -1 && bcIds[direction] != elem.
m_bndCndId) {
3465 ss <<
"Bad cell: " << cellId <<
". Cell has two different boundary condition ids (" << bcIds[direction] <<
", "
3466 << elem.
m_bndCndId <<
") in direction " << direction <<
". This feature is not yet implemented." << endl;
3502 for(
MInt i = 0; i < nDim; i++) {
3503 if(!
approx(normal[i], 0.0, MFloatEps)) {
3504 MInt direction = (normal[i] < 0) ? 2 * i : 2 * i + 1;
3507 direction = 2 * (direction / 2) + 1 - (direction % 2);
3512 array<MFloat, nDim> nghbrCellCoord;
3513 for(
MInt d = 0; d < nDim; ++d) {
3514 nghbrCellCoord[d] = grid().tree().coordinate(cellId, d);
3516 nghbrCellCoord[i] += (2 * (direction % 2) - 1) * grid().cellLengthAtCell(cellId);
3517 const MBool nghbrIsInside = geometry().pointIsInside2(&nghbrCellCoord[0]);
3518 if((!nghbrIsInside && grid().tree().level(cellId) != getLevelOfDirectNeighborCell(cellId, direction))
3519 || (grid().tree().level(cellId) < getLevelOfDirectNeighborCell(cellId, direction))) {
3521 ss <<
"Bad cell: " << cellId <<
". There is h-refinement being used along a non axis-aligned surface in "
3522 <<
"direction " << direction <<
". This feature is not yet implemented." << endl;
3528 if(bcIds[direction] > -1 && bcIds[direction] != elem.
m_bndCndId) {
3530 ss <<
"Bad cell: " << cellId <<
". Cell has two different boundary condition ids (" << bcIds[direction]
3531 <<
", " << elem.
m_bndCndId <<
") in direction " << direction
3532 <<
". This feature is not yet implemented." << endl;
3543 std::vector<CutCandidate<nDim>> cutCandidates(1);
3545 cutCandidates[0].cellId = cellId;
3547 m_geometryIntersection->computeCutPointsFromSTL(cutCandidates);
3549 const MInt noCutPoints = cutCandidates[0].noCutPoints;
3552 vector<MFloat> targetPoints(noCutPoints * nDim);
3553 for(
MInt cutPoint = 0; cutPoint < noCutPoints; cutPoint++) {
3554 for(
MInt i = 0; i < nDim; i++) {
3555 targetPoints[(cutPoint * nDim) + i] = cutCandidates[0].cutPoints[cutPoint][i];
3560 elem.
calcNormal(&targetPoints[0], &normal[0]);
3564 MFloat maxNormal = fabs(normal[0]);
3565 for(
MInt d = 0; d < nDim; d++) {
3566 if(fabs(normal[d]) > maxNormal) {
3568 maxNormal = fabs(normal[d]);
3575 array<MFloat, nDim> centroid;
3577 const MFloat pos = grid().tree().coordinate(cellId, orientation) + 0.5 * grid().cellLengthAtCell(cellId);
3578 const MFloat neg = grid().tree().coordinate(cellId, orientation) - 0.5 * grid().cellLengthAtCell(cellId);
3579 MInt direction = (std::fabs(neg - centroid[orientation]) < std::fabs(pos - centroid[orientation]))
3581 : 2 * orientation + 1;
3587 if((!hasElement && !grid().tree().hasNeighbor(cellId, direction))
3588 || (hasElement && grid().tree().hasNeighbor(cellId, direction))) {
3589 direction = 2 * (direction / 2) + 1 - (direction % 2);
3593 array<MFloat, nDim> nghbrCellCoord;
3594 for(
MInt d = 0; d < nDim; ++d) {
3595 nghbrCellCoord[d] = grid().tree().coordinate(cellId, d);
3597 nghbrCellCoord[orientation] += (2 * (direction % 2) - 1) * grid().cellLengthAtCell(cellId);
3598 const MBool nghbrIsInside = geometry().pointIsInside2(&nghbrCellCoord[0]);
3599 if((!nghbrIsInside && grid().tree().level(cellId) != getLevelOfDirectNeighborCell(cellId, direction))
3600 || (grid().tree().level(cellId) < getLevelOfDirectNeighborCell(cellId, direction))) {
3602 ss <<
"Bad cell: " << cellId <<
". There is h-refinement being used along a non axis-aligned surface "
3603 <<
"in direction " << direction <<
". This feature is not yet implemented." << endl;
3608 if(bcIds[direction] > -1 && bcIds[direction] != elem.
m_bndCndId) {
3610 ss <<
"Bad cell: " << cellId <<
". Cell has two different boundary condition ids (" << bcIds[direction]
3611 <<
", " << elem.
m_bndCndId <<
") in direction " << direction <<
". This feature is not yet implemented."
3643template <MInt nDim,
class SysEqn>
3648 const MInt cellId = m_elements.cellId(elementId);
3651 vector<MInt> nghbrIds;
3652 const MInt parentNghbrId = grid().tree().neighbor(cellId, dir);
3656 static constexpr MInt dirNghbrs[6][4] = {
3665 MInt noNonHaloChilds = 0;
3667 for(
MInt i = 0; i < 2 * (nDim - 1); i++) {
3668 if(grid().tree().hasChild(parentNghbrId, dirNghbrs[dir][i])) {
3670 const MInt childId = grid().tree().child(parentNghbrId, dirNghbrs[dir][i]);
3671 if(!grid().tree().hasProperty(childId, Cell::IsHalo)) {
3672 ASSERT(grid().tree().hasProperty(parentNghbrId, Cell::IsPartLvlAncestor),
3673 "this case is only valid with a neighboring partition level ancestor cell");
3677 nghbrIds.push_back(grid().tree().child(parentNghbrId, dirNghbrs[dir][i]));
3682 if(noNonHaloChilds == 2 * (nDim - 1)) {
3686 if(nghbrIds.empty()) {
3687 stringstream errorMessage;
3688 errorMessage <<
"Error: child cells of neighboring halo cell not found. "
3689 "Try setting 'newCreateWindowCellsF = 0'."
3691 TERMM(1, errorMessage.str());
3695 MInt firstSurfaceId = -1;
3696 for(
const auto& nghbrId : nghbrIds) {
3698 const MInt srfcId = createSurface(elementId, dir, nghbrId);
3701 if(firstSurfaceId == -1) {
3702 firstSurfaceId = srfcId;
3706 for(
MInt hElementId = 0; hElementId < m_helements.size(); hElementId++) {
3708 if(m_helements.elementId(hElementId) == elementId) {
3710 for(
MInt hSurfId = 0; hSurfId < 2 * (nDim - 1); hSurfId++) {
3711 if(m_helements.hrefSurfaceIds(hElementId, dir, hSurfId) == -1) {
3712 m_helements.hrefSurfaceIds(hElementId, dir, hSurfId) = srfcId;
3721 m_surfaces.fineCellId(srfcId) = nghbrId;
3727 return firstSurfaceId;
3749template <MInt nDim,
class SysEqn>
3754 const MInt cellId = m_elements.cellId(elementId);
3755 const MInt parentId = grid().tree().parent(cellId);
3760 nghbrId = grid().tree().neighbor(cellId, dir);
3761 if(nghbrId == -1 && parentId > -1 && grid().tree().hasNeighbor(parentId, dir)) {
3762 nghbrId = grid().tree().neighbor(parentId, dir);
3765 const MInt nghbrElementId = m_elements.getElementByCellId(nghbrId);
3768 const MInt srfcId = m_surfaces.size();
3769 m_surfaces.append();
3772 const MInt oppositeDir = 2 * (dir / 2) + 1 - (dir % 2);
3775 const MInt orientation = dir / 2;
3776 m_surfaces.orientation(srfcId) = orientation;
3780 const MInt nghbrSideId = dir % 2;
3781 const MInt elementSideId = (nghbrSideId + 1) % 2;
3782 m_surfaces.nghbrElementIds(srfcId, nghbrSideId) = nghbrElementId;
3783 m_surfaces.nghbrElementIds(srfcId, elementSideId) = elementId;
3786 m_elements.surfaceIds(elementId, dir) = srfcId;
3789 if(nghbrElementId > -1 && m_elements.surfaceIds(nghbrElementId, oppositeDir) == -1) {
3790 m_elements.surfaceIds(nghbrElementId, oppositeDir) = srfcId;
3794 const MInt level = getLevelByElementId(elementId);
3796 if(nghbrElementId > -1) {
3797 nghbrLvl = getLevelByElementId(nghbrElementId);
3801 m_surfaces.fineCellId(srfcId) = -1;
3804 if(nghbrLvl != -1 && level > nghbrLvl) {
3806 for(
MInt hElementId = 0; hElementId < m_helements.size(); hElementId++) {
3808 if(m_helements.elementId(hElementId) == nghbrElementId) {
3809 m_surfaces.fineCellId(srfcId) = cellId;
3811 for(
MInt hSurfId = 0; hSurfId < 2 * (nDim - 1); hSurfId++) {
3812 if(m_helements.hrefSurfaceIds(hElementId, oppositeDir, hSurfId) == -1) {
3813 m_helements.hrefSurfaceIds(hElementId, oppositeDir, hSurfId) = srfcId;
3824 if(nghbrElementId > -1) {
3825 m_surfaces.polyDeg(srfcId) = max(m_elements.polyDeg(elementId), m_elements.polyDeg(nghbrElementId));
3826 m_surfaces.noNodes1D(srfcId) = max(m_elements.noNodes1D(elementId), m_elements.noNodes1D(nghbrElementId));
3828 m_surfaces.polyDeg(srfcId) = m_elements.polyDeg(elementId);
3829 m_surfaces.noNodes1D(srfcId) = m_elements.noNodes1D(elementId);
3835 nghbrLvl = grid().tree().level(nghbrId);
3839 if(level > nghbrLvl) {
3840 const MFloat cellLength = grid().cellLengthAtCell(cellId);
3841 m_surfaces.coords(srfcId, orientation) =
3842 grid().tree().coordinate(cellId, orientation) + (
static_cast<MFloat>(nghbrSideId) - F1B2) * cellLength;
3843 for(
MInt i = 0; i < nDim; i++) {
3844 if(i != orientation) {
3845 m_surfaces.coords(srfcId, i) = grid().tree().coordinate(cellId, i);
3849 const MFloat cellLength = grid().cellLengthAtCell(nghbrId);
3850 m_surfaces.coords(srfcId, orientation) =
3851 grid().tree().coordinate(nghbrId, orientation) + (
static_cast<MFloat>(1 - nghbrSideId) - F1B2) * cellLength;
3852 for(
MInt i = 0; i < nDim; i++) {
3853 if(i != orientation) {
3854 m_surfaces.coords(srfcId, i) = grid().tree().coordinate(nghbrId, i);
3860 const MInt elementIdL = m_surfaces.nghbrElementIds(srfcId, 0);
3861 const MInt elementIdR = m_surfaces.nghbrElementIds(srfcId, 1);
3864 m_surfaces.internalSideId(srfcId) = -1;
3865 if(elementIdL == -1) {
3866 m_surfaces.internalSideId(srfcId) = 1;
3868 if(elementIdR == -1) {
3869 m_surfaces.internalSideId(srfcId) = 0;
3873 calcSurfaceNodeCoordinates(srfcId);
3876 const MLong globalCellId = grid().tree().globalId(m_elements.cellId(elementId));
3877 if(nghbrId == -1 || nghbrLvl == -1) {
3879 m_surfaces.globalId(srfcId) = 2 * nDim * globalCellId + dir;
3883 const MLong globalNghbrId = grid().tree().globalId(nghbrId);
3884 if(level > nghbrLvl || (globalCellId < globalNghbrId && level == nghbrLvl)) {
3885 m_surfaces.globalId(srfcId) = 2 * nDim * globalCellId + dir;
3887 m_surfaces.globalId(srfcId) = 2 * nDim * globalNghbrId + oppositeDir;
3904template <MInt nDim,
class SysEqn>
3906 const MInt internalSide = (m_surfaces.internalSideId(srfcId) <= 0) ? 0 : 1;
3907 const MInt elementId = m_surfaces.nghbrElementIds(srfcId, internalSide);
3908 const MInt cellId = m_elements.cellId(elementId);
3911 const MInt polyDeg = m_surfaces.polyDeg(srfcId);
3912 const MInt noNodes1D = m_surfaces.noNodes1D(srfcId);
3913 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
3914 const MFloat length = grid().cellLengthAtCell(cellId);
3915 MFloatTensor nodeCoordinates(&m_surfaces.nodeCoords(srfcId), noNodes1D, noNodes1D3, nDim);
3920 switch(m_surfaces.orientation(srfcId)) {
3922 for(
MInt j = 0; j < noNodes1D; j++) {
3923 for(
MInt k = 0; k < noNodes1D3; k++) {
3924 nodeCoordinates(j, k, 0) = m_surfaces.coords(srfcId, 0);
3925 nodeCoordinates(j, k, 1) =
3926 m_surfaces.coords(srfcId, 1) + F1B2 * length * m_interpolation[polyDeg][noNodes1D].m_nodes[j];
3927 IF_CONSTEXPR(nDim == 3) {
3928 nodeCoordinates(j, k, 2) =
3929 m_surfaces.coords(srfcId, 2) + F1B2 * length * m_interpolation[polyDeg][noNodes1D].m_nodes[k];
3936 for(
MInt i = 0; i < noNodes1D; i++) {
3937 for(
MInt k = 0; k < noNodes1D3; k++) {
3938 nodeCoordinates(i, k, 0) =
3939 m_surfaces.coords(srfcId, 0) + F1B2 * length * m_interpolation[polyDeg][noNodes1D].m_nodes[i];
3940 nodeCoordinates(i, k, 1) = m_surfaces.coords(srfcId, 1);
3941 IF_CONSTEXPR(nDim == 3) {
3942 nodeCoordinates(i, k, 2) =
3943 m_surfaces.coords(srfcId, 2) + F1B2 * length * m_interpolation[polyDeg][noNodes1D].m_nodes[k];
3950 for(
MInt i = 0; i < noNodes1D; i++) {
3951 for(
MInt j = 0; j < noNodes1D3; j++) {
3952 nodeCoordinates(i, j, 0) =
3953 m_surfaces.coords(srfcId, 0) + F1B2 * length * m_interpolation[polyDeg][noNodes1D].m_nodes[i];
3954 nodeCoordinates(i, j, 1) =
3955 m_surfaces.coords(srfcId, 1) + F1B2 * length * m_interpolation[polyDeg][noNodes1D].m_nodes[j];
3956 IF_CONSTEXPR(nDim == 3) { nodeCoordinates(i, j, 2) = m_surfaces.coords(srfcId, 2); }
3962 mTerm(1, AT_,
"Bad orientation");
3971template <MInt nDim,
class SysEqn>
3979 const MBool neighbors = m_isInitMpiExchange ? m_noExchangeNghbrDomains > 0 : grid().noNeighborDomains() > 0;
3983 return (noDomains() > 1 || neighbors);
3995template <MInt nDim,
class SysEqn>
3999 m_log <<
"Initialize MPI exchange... ";
4002 map<MInt, MInt> exchangeNghbrDomainMap;
4007 m_noExchangeNghbrDomains = 0;
4010 m_exchangeNghbrDomains.resize(grid().noNeighborDomains());
4011 m_mpiSurfaces.resize(grid().noNeighborDomains());
4015 const MInt begin = m_mpiSurfacesOffset;
4016 const MInt end = m_mpiSurfacesOffset + m_noMpiSurfaces;
4017 for(
MInt srfcId = begin; srfcId < end; srfcId++) {
4018 const MInt internalSideId = m_surfaces.internalSideId(srfcId);
4019 const MInt elementId = m_surfaces.nghbrElementIds(srfcId, internalSideId);
4020 const MInt internalCellId = m_elements.cellId(elementId);
4021 const MInt nghbrCellDir = 2 * m_surfaces.orientation(srfcId) + 1 - internalSideId;
4024 MInt haloCellId = grid().tree().neighbor(internalCellId, nghbrCellDir);
4025 if(haloCellId == -1) {
4027 haloCellId = grid().tree().neighbor(grid().tree().parent(internalCellId), nghbrCellDir);
4031 if(grid().tree().hasChildren(haloCellId)) {
4032 haloCellId = m_surfaces.fineCellId(srfcId);
4036 MInt exchangeNghbrDomainId = -1;
4037 for(
MInt i = 0; i < noDomains() + 1; i++) {
4038 if(grid().domainOffset(i) > grid().tree().globalId(haloCellId)) {
4039 exchangeNghbrDomainId = i - 1;
4043 ASSERT(exchangeNghbrDomainId > -1,
"Could not find domain id!");
4046 if(exchangeNghbrDomainMap.count(exchangeNghbrDomainId) == 0) {
4047 m_exchangeNghbrDomains[m_noExchangeNghbrDomains] = exchangeNghbrDomainId;
4048 exchangeNghbrDomainMap[exchangeNghbrDomainId] = m_noExchangeNghbrDomains;
4049 m_noExchangeNghbrDomains++;
4053 m_mpiSurfaces[exchangeNghbrDomainMap[exchangeNghbrDomainId]].push_back(srfcId);
4057 m_exchangeNghbrDomains.resize(m_noExchangeNghbrDomains);
4058 m_mpiSurfaces.resize(m_noExchangeNghbrDomains);
4061 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4062 sort(m_mpiSurfaces[i].begin(), m_mpiSurfaces[i].end(),
4063 [
this](
const MInt a,
const MInt b) {
return (m_surfaces.globalId(
a) < m_surfaces.globalId(
b)); });
4067 m_mpiRecvSurfaces.resize(m_noExchangeNghbrDomains);
4068 m_mpiRecvSurfaces = m_mpiSurfaces;
4073 if(grid().periodicCartesianDir(0) || grid().periodicCartesianDir(1)
4074 || (nDim == 3 && grid().periodicCartesianDir(2))) {
4075 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4077 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size() - 1; j++) {
4078 if(m_surfaces.globalId(m_mpiSurfaces[i][j]) == m_surfaces.globalId(m_mpiSurfaces[i][j + 1])) {
4079 swap(m_mpiRecvSurfaces[i][j], m_mpiRecvSurfaces[i][j + 1]);
4111#ifdef DG_USE_MPI_BUFFERS
4113 m_sendBuffers.resize(m_noExchangeNghbrDomains);
4114 m_recvBuffers.resize(m_noExchangeNghbrDomains);
4115 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4117 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
4119 const MInt noNodes1D = m_maxNoNodes1D;
4120 const MInt noNodesXD =
ipow(noNodes1D, nDim - 1);
4121 const MInt dataBlockSize = noNodesXD * SysEqn::noVars();
4122 size += dataBlockSize;
4125 m_sendBuffers[i].resize(size);
4126 m_recvBuffers[i].resize(size);
4129#ifdef DG_USE_MPI_DERIVED_TYPES
4130#error Does not work anymore - need to fix for p-refinement!
4132 m_sendTypes.resize(m_noExchangeNghbrDomains);
4133 m_recvTypes.resize(m_noExchangeNghbrDomains);
4134 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4138 const MInt count = m_mpiSurfaces[i].size();
4144 m_sendTypes[i] = MPI_DATATYPE_NULL;
4145 m_recvTypes[i] = MPI_DATATYPE_NULL;
4149 vector<MInt> lengths(count);
4150 vector<MPI_Aint> displacementsSend(count);
4151 vector<MPI_Aint> displacementsRecv(count);
4154 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
4155 const MInt srfcId = m_mpiSurfaces[i][j];
4157 const MInt noNodes1D = m_surfaces.noNodes1D(srfcId);
4158 const MInt noNodesXD =
ipow(noNodes1D, nDim - 1);
4159 const MInt dataBlockSize = noNodesXD * SysEqn::noVars();
4160 lengths[j] = dataBlockSize;
4163 const MInt internalSideId = m_surfaces.internalSideId(srfcId);
4164 MPI_Get_address(m_surfaces.variables(srfcId, internalSideId), &displacementsSend[j], AT_);
4165 MPI_Get_address(m_surfaces.variables(srfcId, 1 - internalSideId), &displacementsRecv[j], AT_);
4182 m_sendRequests.resize(m_noExchangeNghbrDomains);
4183 m_recvRequests.resize(m_noExchangeNghbrDomains);
4184 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4185#ifdef DG_USE_MPI_BUFFERS
4187 m_exchangeNghbrDomains[i], domainId(), mpiComm(), &m_sendRequests[i], AT_,
"m_sendBuffers[i][0]");
4189 m_exchangeNghbrDomains[i], m_exchangeNghbrDomains[i], mpiComm(), &m_recvRequests[i], AT_,
4190 "m_recvBuffers[i][0]");
4192#ifdef DG_USE_MPI_DERIVED_TYPES
4193 MPI_Send_init(MPI_BOTTOM, 1, m_sendTypes[i], m_exchangeNghbrDomains[i], domainId(), mpiComm(), &m_sendRequests[i],
4195 MPI_Recv_init(MPI_BOTTOM, 1, m_recvTypes[i], m_exchangeNghbrDomains[i], m_exchangeNghbrDomains[i], mpiComm(),
4196 &m_recvRequests[i], AT_,
"MPI_BOTTOM");
4200 m_isInitMpiExchange =
true;
4201 m_log <<
"done" << endl;
4214template <MInt nDim,
class SysEqn>
4218 m_log <<
"Initializing simulation statistics... ";
4223 m_statLocalMinLevel = grid().minLevel();
4224 m_statLocalMaxLevel = grid().maxLevel();
4227 m_statLocalNoCells = grid().noCells();
4228 m_statLocalNoInternalCells = grid().noInternalCells();
4229 m_statLocalNoHaloCells = grid().noCells() - grid().noInternalCells();
4230 m_statLocalMaxNoCells = grid().raw().treeb().capacity();
4233 m_statLocalNoElements = m_elements.size();
4234 m_statLocalNoHElements = m_helements.size();
4237 m_statLocalNoSurfaces = m_surfaces.size();
4238 m_statLocalNoBoundarySurfaces = m_noBoundarySurfaces;
4239 m_statLocalNoInnerSurfaces = m_noInnerSurfaces;
4240 m_statLocalNoMpiSurfaces = m_noMpiSurfaces;
4241 m_statLocalMaxNoSurfaces = m_maxNoSurfaces;
4244 m_statLocalNoActiveCells = m_elements.size();
4245 m_statLocalNoActiveDOFs =
4246 accumulate(&m_elements.noNodes1D(0), &m_elements.noNodes1D(0) + m_elements.size(), 0,
4247 [](
const MInt result,
const MInt noNodes1D) { return result + ipow(noNodes1D, nDim); });
4250 const MInt noPolyDegs = m_maxPolyDeg - m_minPolyDeg + 1;
4251 if(noPolyDegs > 1) {
4252 m_statLocalNoActiveDOFsPolyDeg.assign(noPolyDegs, 0);
4253 const MInt noElements = m_elements.size();
4254 for(
MInt elementId = 0; elementId < noElements; elementId++) {
4255 const MInt polyDeg = m_elements.polyDeg(elementId);
4256 const MInt noDOFs = m_elements.noNodesXD(elementId);
4257 m_statLocalNoActiveDOFsPolyDeg[polyDeg - m_minPolyDeg] += noDOFs;
4263 m_statLocalMinPolyDeg = *min_element(&m_elements.polyDeg(0), &m_elements.polyDeg(0) + m_elements.size());
4264 m_statLocalMaxPolyDeg = *max_element(&m_elements.polyDeg(0), &m_elements.polyDeg(0) + m_elements.size());
4267 m_statLocalPolyDegSum = accumulate(&m_elements.polyDeg(0), &m_elements.polyDeg(0) + m_elements.size(), 0);
4270 m_statLocalNoHrefSurfs = count_if(&m_surfaces.fineCellId(0), &m_surfaces.fineCellId(0) + m_surfaces.size(),
4271 [](
const MInt id) { return id != -1; });
4272 m_statLocalNoPrefSurfs = 0;
4275 for(
MInt srfcId = m_innerSurfacesOffset; srfcId < m_innerSurfacesOffset + m_noInnerSurfaces; srfcId++) {
4276 const MInt elementIdL = m_surfaces.nghbrElementIds(srfcId, 0);
4277 const MInt elementIdR = m_surfaces.nghbrElementIds(srfcId, 1);
4278 if(m_elements.polyDeg(elementIdL) != m_elements.polyDeg(elementIdR)) {
4279 m_statLocalNoPrefSurfs++;
4288 buffer[0] = m_statLocalMinLevel;
4289 buffer[1] = m_statLocalMinPolyDeg;
4291 buffer[2] = m_statLocalMaxLevel;
4292 buffer[3] = m_statLocalMaxPolyDeg;
4294 buffer[4] = m_statLocalNoCells;
4295 buffer[5] = m_statLocalNoInternalCells;
4296 buffer[6] = m_statLocalNoHaloCells;
4297 buffer[7] = m_statLocalMaxNoCells;
4298 buffer[8] = m_statLocalNoElements;
4299 buffer[9] = m_statLocalNoSurfaces;
4300 buffer[10] = m_statLocalNoBoundarySurfaces;
4301 buffer[11] = m_statLocalNoInnerSurfaces;
4302 buffer[12] = m_statLocalNoMpiSurfaces;
4303 buffer[13] = m_statLocalMaxNoSurfaces;
4304 buffer[14] = m_statLocalNoActiveCells;
4305 buffer[15] = m_statLocalNoActiveDOFs;
4306 buffer[16] = m_statLocalPolyDegSum;
4307 buffer[17] = m_statLocalNoHrefSurfs;
4308 buffer[18] = m_statLocalNoPrefSurfs;
4309 buffer[19] = m_statLocalNoHElements;
4310 buffer[20] = m_noCutOffBoundarySurfaces[0];
4311 buffer[21] = m_noCutOffBoundarySurfaces[1];
4312 buffer[22] = m_noCutOffBoundarySurfaces[2];
4313 buffer[23] = m_noCutOffBoundarySurfaces[3];
4314 IF_CONSTEXPR(nDim == 3) {
4315 buffer[24] = m_noCutOffBoundarySurfaces[4];
4316 buffer[25] = m_noCutOffBoundarySurfaces[5];
4323 RECORD_TIMER_START(m_timers[Timers::Accumulated]);
4324 RECORD_TIMER_START(m_timers[Timers::MPI]);
4325 RECORD_TIMER_START(m_timers[Timers::MPIComm]);
4329 MPI_Allreduce(MPI_IN_PLACE, &buffer[0], 2, MPI_INT, MPI_MIN, mpiComm(), AT_,
"MPI_IN_PLACE",
"buffer[0]");
4331 MPI_Allreduce(MPI_IN_PLACE, &buffer[2], 2, MPI_INT, MPI_MAX, mpiComm(), AT_,
"MPI_IN_PLACE",
"buffer[2]");
4334 "MPI_IN_PLACE",
"buffer[4]");
4336 RECORD_TIMER_STOP(m_timers[Timers::MPIComm]);
4337 RECORD_TIMER_STOP(m_timers[Timers::MPI]);
4338 RECORD_TIMER_STOP(m_timers[Timers::Accumulated]);
4341 m_statGlobalMinLevel = buffer[0];
4342 m_statGlobalMinPolyDeg = buffer[1];
4343 m_statGlobalMaxLevel = buffer[2];
4344 m_statGlobalMaxPolyDeg = buffer[3];
4345 m_statGlobalNoCells = buffer[4];
4346 m_statGlobalNoInternalCells = buffer[5];
4347 m_statGlobalNoHaloCells = buffer[6];
4348 m_statGlobalMaxNoCells = buffer[7];
4349 m_statGlobalNoElements = buffer[8];
4350 m_statGlobalNoSurfaces = buffer[9];
4351 m_statGlobalNoBoundarySurfaces = buffer[10];
4352 m_statGlobalNoInnerSurfaces = buffer[11];
4353 m_statGlobalNoMpiSurfaces = buffer[12];
4354 m_statGlobalMaxNoSurfaces = buffer[13];
4355 m_statGlobalNoActiveCells = buffer[14];
4356 m_statGlobalNoActiveDOFs = buffer[15];
4357 m_statGlobalPolyDegSum = buffer[16];
4358 m_statGlobalNoHrefSurfs = buffer[17];
4359 m_statGlobalNoPrefSurfs = buffer[18];
4360 m_statGlobalNoHElements = buffer[19];
4361 m_statGlobalNoCutOffBoundarySurfaces[0] = buffer[20];
4362 m_statGlobalNoCutOffBoundarySurfaces[1] = buffer[21];
4363 m_statGlobalNoCutOffBoundarySurfaces[2] = buffer[22];
4364 m_statGlobalNoCutOffBoundarySurfaces[3] = buffer[23];
4365 m_statGlobalNoCutOffBoundarySurfaces[4] = buffer[24];
4366 m_statGlobalNoCutOffBoundarySurfaces[5] = buffer[25];
4368 m_log <<
"done" << endl;
4383template <MInt nDim,
class SysEqn>
4389 RECORD_TIMER_STOP(m_timers[Timers::Run]);
4393 RECORD_TIMER_START(m_timers[Timers::CleanUp]);
4396 if(!m_isInitSolver) {
4397 TERMM(1,
"Solver was not initialized.");
4401 if(hasMpiExchange()) {
4402 finalizeMpiExchange();
4405 RECORD_TIMER_STOP(m_timers[Timers::CleanUp]);
4406 RECORD_TIMER_STOP(m_timers[Timers::Run]);
4417template <MInt nDim,
class SysEqn>
4423 RECORD_TIMER_START(m_timers[Timers::MainLoop]);
4424 RECORD_TIMER_START(m_timers[Timers::RungeKuttaStep]);
4425 RECORD_TIMER_START(m_timers[Timers::TimeDeriv]);
4426 finishMpiSurfaceExchange();
4427 RECORD_TIMER_STOP(m_timers[Timers::TimeDeriv]);
4428 RECORD_TIMER_STOP(m_timers[Timers::RungeKuttaStep]);
4429 RECORD_TIMER_STOP(m_timers[Timers::MainLoop]);
4432 MPI_Waitall(m_noExchangeNghbrDomains, &m_sendRequests[0], MPI_STATUSES_IGNORE, AT_);
4433 m_mpiSendRequestsOpen =
false;
4437 if(m_mpiRecvRequestsOpen) {
4438 std::vector<MBool> waitForCancel(m_noExchangeNghbrDomains,
false);
4439 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4440 if(m_recvRequests[i] != MPI_REQUEST_NULL) {
4442 waitForCancel[i] =
true;
4446 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4447 if(waitForCancel[i]) {
4448 MPI_Wait(&m_recvRequests[i], MPI_STATUS_IGNORE, AT_);
4451 m_mpiRecvRequestsOpen =
false;
4463template <MInt nDim,
class SysEqn>
4467 if(m_noExchangeNghbrDomains > 0) {
4469 if(m_mpiSendRequestsOpen) {
4470 MPI_Waitall(m_noExchangeNghbrDomains, &m_sendRequests[0], MPI_STATUSES_IGNORE, AT_);
4474 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4475 if(m_sendRequests[i] != MPI_REQUEST_NULL) {
4482 if(m_recvRequests[i] != MPI_REQUEST_NULL) {
4484 if(m_mpiRecvRequestsOpen) {
4492 m_mpiSendRequestsOpen =
false;
4493 m_mpiRecvRequestsOpen =
false;
4495#ifdef DG_USE_MPI_DERIVED_TYPES
4497 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4511template <MInt nDim,
class SysEqn>
4514 RECORD_TIMER_START(m_timers[Timers::InitialCondition]);
4516 m_log <<
"Applying initial conditions... ";
4518 const MFloat time = m_startTime;
4525 const MInt noElements = m_elements.size();
4526 for(
MInt elementId = 0; elementId < noElements; elementId++) {
4527 const MInt noNodes1D = m_elements.noNodes1D(elementId);
4528 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
4531 MFloatTensor nodeVars(&m_elements.nodeVars(elementId), noNodes1D, noNodes1D, noNodes1D3,
4532 max(SysEqn::noNodeVars(), 1));
4533 MFloatTensor u(&m_elements.variables(elementId), noNodes1D, noNodes1D, noNodes1D3, SysEqn::noVars());
4534 MFloatTensor x(&m_elements.nodeCoordinates(elementId), noNodes1D, noNodes1D, noNodes1D3, nDim);
4536 for(
MInt i = 0; i < noNodes1D; i++) {
4537 for(
MInt j = 0; j < noNodes1D; j++) {
4538 for(
MInt k = 0; k < noNodes1D3; k++) {
4540 m_sysEqn.calcInitialCondition(time, &x(i, j, k, 0), &nodeVars(i, j, k, 0), &u(i, j, k, 0));
4546 m_log <<
"done" << endl;
4548 RECORD_TIMER_STOP(m_timers[Timers::InitialCondition]);
4564template <MInt nDim,
class SysEqn>
4568 m_log <<
"Update node variable data on surfaces... ";
4570 const MInt*
const polyDegs = &m_elements.polyDeg(0);
4571 const MInt*
const noNodes = &m_elements.noNodes1D(0);
4572 const MInt*
const surfaceIds = &m_elements.surfaceIds(0, 0);
4573 const MInt noElements = m_elements.size();
4583 using namespace dg::interpolation;
4587#pragma omp parallel for
4590 for(
MInt elementId = 0; elementId < noElements; elementId++) {
4591 const MInt surfaceIdOffset = elementId * 2 * nDim;
4592 const MInt polyDeg = polyDegs[elementId];
4593 const MInt noNodes1D = noNodes[elementId];
4597 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
4598 const MInt srfcId = surfaceIds[surfaceIdOffset + dir];
4599 const MInt side = 1 - dir % 2;
4601 MFloat* src = &m_elements.nodeVars(elementId);
4602 MFloat* dest = &m_surfaces.nodeVars(srfcId, side);
4603 prolongToFaceGauss<nDim, SysEqn::noNodeVars()>(src, dir, noNodes1D, &interp.
m_LFace[0][0],
4609#pragma omp parallel for
4612 for(
MInt elementId = 0; elementId < noElements; elementId++) {
4613 const MInt surfaceIdOffset = elementId * 2 * nDim;
4614 const MInt noNodes1D = noNodes[elementId];
4617 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
4618 const MInt srfcId = surfaceIds[surfaceIdOffset + dir];
4619 const MInt side = 1 - dir % 2;
4621 MFloat* src = &m_elements.nodeVars(elementId);
4622 MFloat* dest = &m_surfaces.nodeVars(srfcId, side);
4623 prolongToFaceGaussLobatto<nDim, SysEqn::noNodeVars()>(src, dir, noNodes1D, dest);
4637 const MInt noVars = SysEqn::noNodeVars();
4638 const MInt maxNoNodes1D = m_maxNoNodes1D;
4639 const MInt maxNoNodes1D3 = (nDim == 3) ? maxNoNodes1D : 1;
4640 const MInt noHElements = m_helements.size();
4641 const MInt noDirs = 2 * nDim;
4642 const MInt noSurfs = 2 * (nDim - 1);
4646 if(noHElements > 0) {
4648#pragma omp parallel for
4650 for(
MInt hElementId = 0; hElementId < noHElements; hElementId++) {
4651 const MInt coarseElementId = m_helements.elementId(hElementId);
4652 for(
MInt dir = 0; dir < noDirs; dir++) {
4653 const MInt coarseSrfcId = m_elements.surfaceIds(coarseElementId, dir);
4654 for(
MInt pos = 0; pos < noSurfs; pos++) {
4655 const MInt hSrfcId = m_helements.hrefSurfaceIds(hElementId, dir, pos);
4656 const MInt side = 1 - dir % 2;
4664 if(coarseSrfcId != hSrfcId) {
4665 const MInt noNodes1D = noNodes[coarseElementId];
4666 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
4667 const MInt noNodesXD = noNodes1D * noNodes1D3;
4670 copy_n(&m_surfaces.nodeVars(coarseSrfcId, side),
4671 noNodesXD * SysEqn::noNodeVars(),
4672 &m_surfaces.nodeVars(hSrfcId, side));
4680 MFloatTensor projected(maxNoNodes1D, maxNoNodes1D3, noVars);
4683 if(noHElements > 0) {
4685#pragma omp parallel for firstprivate(projected)
4687 for(
MInt hElementId = 0; hElementId < noHElements; hElementId++) {
4688 for(
MInt dir = 0; dir < noDirs; dir++) {
4689 for(
MInt pos = 0; pos < noSurfs; pos++) {
4690 const MInt hSrfcId = m_helements.hrefSurfaceIds(hElementId, dir, pos);
4691 const MInt side = 1 - dir % 2;
4698 const MInt surfaceNoNodes1D = m_surfaces.noNodes1D(hSrfcId);
4699 const MInt surfaceNoNodes1D3 = (nDim == 3) ? surfaceNoNodes1D : 1;
4700 const MInt size = surfaceNoNodes1D * surfaceNoNodes1D3 * noVars;
4703 calcMortarProjection<dg::mortar::forward, SysEqn::noNodeVars()>(
4704 hSrfcId, dir, &m_surfaces.nodeVars(hSrfcId, side), &projected[0], m_elements, m_surfaces);
4707 copy_n(&projected[0], size, &m_surfaces.nodeVars(hSrfcId, side));
4715#pragma omp parallel for firstprivate(projected)
4717 for(
MInt elementId = 0; elementId < noElements; elementId++) {
4718 const MInt surfaceIdOffset = elementId * 2 * nDim;
4720 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
4722 const MInt srfcId = surfaceIds[surfaceIdOffset + dir];
4725 if(srfcId < m_innerSurfacesOffset) {
4729 const MInt surfacePolyDeg = m_surfaces.polyDeg(srfcId);
4730 const MInt elementPolyDeg = m_elements.polyDeg(elementId);
4731 const MInt side = 1 - dir % 2;
4732 const MInt surfaceNoNodes1D = m_surfaces.noNodes1D(srfcId);
4733 const MInt surfaceNoNodes1D3 = (nDim == 3) ? surfaceNoNodes1D : 1;
4734 const MInt size = surfaceNoNodes1D * surfaceNoNodes1D3 * noVars;
4737 if(m_surfaces.fineCellId(srfcId) != -1 && m_surfaces.fineCellId(srfcId) != m_elements.cellId(elementId)) {
4742 if(surfacePolyDeg > elementPolyDeg) {
4744 calcMortarProjection<dg::mortar::forward, SysEqn::noNodeVars()>(srfcId, dir, &m_surfaces.nodeVars(srfcId, side),
4745 &projected[0], m_elements, m_surfaces);
4748 copy_n(&projected[0], size, &m_surfaces.nodeVars(srfcId, side));
4758 if(!hasMpiExchange()) {
4768 vector<vector<MFloat>> sendBuffers(m_noExchangeNghbrDomains);
4769 vector<vector<MFloat>> recvBuffers(m_noExchangeNghbrDomains);
4770 const MInt dataBlockSize =
ipow(m_maxNoNodes1D, nDim - 1) * SysEqn::noNodeVars();
4771 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4772 const MInt size = m_mpiSurfaces[i].size() * dataBlockSize;
4773 sendBuffers[i].resize(size);
4774 recvBuffers[i].resize(size);
4777 fill(sendRequests.
begin(), sendRequests.
end(), MPI_REQUEST_NULL);
4779 fill(recvRequests.
begin(), recvRequests.
end(), MPI_REQUEST_NULL);
4782 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4784 m_exchangeNghbrDomains[i], mpiComm(), &recvRequests[i], AT_,
"recvBuffers[i][0]");
4788 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4790 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
4791 const MInt srfcId = m_mpiSurfaces[i][j];
4792 const MInt sideId = m_surfaces.internalSideId(srfcId);
4795 const MFloat* data = &m_surfaces.nodeVars(srfcId, sideId);
4796 copy_n(data, dataBlockSize, &sendBuffers[i][size]);
4797 size += dataBlockSize;
4799 ASSERT(size ==
static_cast<MInt>(sendBuffers[i].size()),
"Data size does not match buffer size.");
4803 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4805 domainId(), mpiComm(), &sendRequests[i], AT_,
"sendBuffers[i][0]");
4809 MPI_Waitall(m_noExchangeNghbrDomains, &recvRequests[0], MPI_STATUSES_IGNORE, AT_);
4812 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
4814 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
4815 const MInt srfcId = m_mpiSurfaces[i][j];
4816 const MInt sideId = 1 - m_surfaces.internalSideId(srfcId);
4819 MFloat* data = &m_surfaces.nodeVars(srfcId, sideId);
4820 std::copy_n(&recvBuffers[i][size], dataBlockSize, data);
4821 size += dataBlockSize;
4823 ASSERT(size ==
static_cast<MInt>(recvBuffers[i].size()),
"Data size does not match buffer size.");
4827 MPI_Waitall(m_noExchangeNghbrDomains, &sendRequests[0], MPI_STATUSES_IGNORE, AT_);
4829 m_log <<
"done" << endl;
4844template <MInt nDim,
class SysEqn>
4847 MBool nodeVarsExtension =
false;
4849 nodeVarsExtension = Context::getSolverProperty<MBool>(
"nodeVarsExtension", m_solverId, AT_, &nodeVarsExtension);
4854 nodeVarsExtension = extensionSetImplicit;
4858 if(!nodeVarsExtension) {
4866 if(noExtendOffsets != noExtendDirs || noExtendLimits != noExtendDirs) {
4867 TERMM(1,
"ERROR: # of specified directions (" + to_string(noExtendDirs) +
"), # of offsets ("
4868 + to_string(noExtendOffsets) +
"), and # of limits (" + to_string(noExtendLimits) +
") do not match!");
4871 std::vector<MInt> meanExtendDirections(noExtendDirs);
4872 std::vector<MFloat> meanExtendOffsets(noExtendDirs);
4873 std::vector<MFloat> meanExtendLimits(noExtendDirs);
4875 for(
MInt i = 0; i < noExtendDirs; i++) {
4876 meanExtendDirections[i] = Context::getSolverProperty<MInt>(
"meanExtendDirections", m_solverId, AT_, i);
4877 meanExtendOffsets[i] = Context::getSolverProperty<MFloat>(
"meanExtendOffsets", m_solverId, AT_, i);
4878 meanExtendLimits[i] = Context::getSolverProperty<MFloat>(
"meanExtendLimits", m_solverId, AT_, i);
4882 for(
MInt i = 0; i < noExtendDirs; i++) {
4883 if(meanExtendDirections[i] >= 0 && meanExtendDirections[i] < 2 * nDim) {
4884 extendNodeVariablesSingleDirection(meanExtendDirections[i], meanExtendOffsets[i], meanExtendLimits[i]);
4903template <MInt nDim,
class SysEqn>
4905 const MFloat extendOffset,
4906 const MFloat extendLimit) {
4907 m_log <<
"Propagating mean values in direction " << extendDir <<
" from offset " << extendOffset <<
"..."
4910 const MInt noElements = m_elements.size();
4911 const MInt*
const noNodes1D = &m_elements.noNodes1D(0);
4912 const MInt extendSide = extendDir % 2;
4913 const MInt extendDimension = (extendDir - extendSide) / 2;
4914 const MInt donorDir = 2 * extendDimension + (1 - extendSide);
4915 const MInt noHElements = m_helements.size();
4916 const MInt noSurfs = 2 * (nDim - 1);
4929 MIntTensor hForwardSurfaceId(m_surfaces.size(), noSurfs);
4930 MIntTensor hReverseSurfaceId(m_surfaces.size(), noSurfs);
4931 hForwardSurfaceId.
set(-1);
4932 hReverseSurfaceId.
set(-1);
4933 for(
MInt hElemId = 0; hElemId < noHElements; hElemId++) {
4934 const MInt elemId = m_helements.elementId(hElemId);
4936 const MInt surfIdForward = m_elements.surfaceIds(elemId, extendDir);
4937 const MInt surfIdReverse = m_elements.surfaceIds(elemId, donorDir);
4938 for(
MInt pos = 0; pos < noSurfs; pos++) {
4939 hForwardSurfaceId(surfIdForward, pos) = m_helements.hrefSurfaceIds(hElemId, extendDir, pos);
4940 const MInt hReverseSurfaceIdFine = m_helements.hrefSurfaceIds(hElemId, donorDir, pos);
4941 hReverseSurfaceId(surfIdReverse, pos) = hReverseSurfaceIdFine;
4943 if(hReverseSurfaceId(surfIdReverse, pos) != surfIdReverse && hReverseSurfaceIdFine > 0) {
4944 hReverseSurfaceId(hReverseSurfaceIdFine, 0) = m_surfaces.
size() + 1;
4949 MIntTensor surfaceHasRecvd(m_surfaces.size());
4950 surfaceHasRecvd.
set(0.0);
4952 std::vector<MBool> elementUpdated(noElements,
false);
4953 std::vector<MBool> elementOutsideLimit(noElements,
false);
4954 std::vector<MInt> mpiSurfaceSendSize(m_surfaces.size(), 0);
4955 for(
MInt elemId = 0; elemId < noElements; elemId++) {
4957 const MFloat pSideCoordinate =
4958 m_surfaces.coords(m_elements.surfaceIds(elemId, extendDimension * 2 + 1), extendDimension);
4959 const MFloat mSideCoordinate =
4960 m_surfaces.coords(m_elements.surfaceIds(elemId, extendDimension * 2), extendDimension);
4961 elementUpdated[elemId] = (pSideCoordinate >= extendOffset && mSideCoordinate <= extendOffset);
4963 if((!extendSide && pSideCoordinate <= extendLimit && mSideCoordinate <= extendLimit)
4964 || (extendSide && pSideCoordinate >= extendLimit && mSideCoordinate >= extendLimit)) {
4965 elementOutsideLimit[elemId] =
true;
4968 if(elementUpdated[elemId]) {
4969 const MInt sId = m_elements.surfaceIds(elemId, extendDir);
4970 if(sId >= m_mpiSurfacesOffset) {
4971 const MInt buffSize =
ipow(noNodes1D[elemId], nDim - 1) * SysEqn::noNodeVars();
4972 mpiSurfaceSendSize[sId] = buffSize;
4973 if(hForwardSurfaceId(sId, 0) > -1) {
4974 for(
MInt pos = 0; pos < noSurfs; pos++) {
4975 const MInt hElemId = m_surfaces.nghbrElementIds(sId, extendSide);
4976 mpiSurfaceSendSize[m_elements.surfaceIds(hElemId, extendDir)] = buffSize;
4983 for(
MInt extendIterator = 0; extendIterator < 4 * noDomains() + 1; extendIterator++) {
4984 if(extendIterator == 4 * noDomains()) {
4986 for(
MInt elemId = 0; elemId < noElements; elemId++) {
4987 if(elementUpdated[elemId] == 0) {
4988 m_log <<
"r" << domainId() <<
" e" << elemId <<
" not updated" << std::endl;
4991 TERMM(1,
"Maximum amount of extension iterations reached.");
4993 for(
MInt elementCounter = 0; elementCounter < noElements; elementCounter++) {
4994 MBool noMoreUpdates =
true;
4995 for(
MInt elemId = 0; elemId < noElements; elemId++) {
4998 const MInt srfcId = m_elements.surfaceIds(elemId, extendDir);
4999 if(srfcId < m_innerSurfacesOffset) {
5002 if(srfcId < m_mpiSurfacesOffset) {
5003 const MInt updateeElementId = m_surfaces.nghbrElementIds(srfcId, extendSide);
5004 if(elementUpdated[elemId] && !(elementUpdated[updateeElementId])
5005 && !(elementOutsideLimit[updateeElementId])) {
5006 updateNodeVariablesSingleElement(updateeElementId, extendDir, hForwardSurfaceId, hReverseSurfaceId,
5007 elementUpdated, mpiSurfaceSendSize, noMoreUpdates);
5008 noMoreUpdates =
false;
5018 MInt extendDomainId;
5019 MPI_Comm_rank(mpiComm(), &extendDomainId);
5021 MInt thisRankIsDone = 1;
5022 for(
MInt a = 0;
a < m_surfaces.size();
a++) {
5023 if(mpiSurfaceSendSize[
a] > 0) {
5027 MInt allRanksAreDone;
5028 MPI_Allreduce(&thisRankIsDone, &allRanksAreDone, 1, MPI_INT, MPI_MIN, mpiComm(), AT_,
"thisRankIsDone",
5030 if(allRanksAreDone) {
5036 std::vector<std::vector<MInt>> nghbrSendSizes(m_noExchangeNghbrDomains);
5037 std::vector<std::vector<MInt>> nghbrRecvSizes(m_noExchangeNghbrDomains);
5039 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
5040 nghbrSendSizes[i].resize(m_mpiSurfaces[i].size());
5041 nghbrRecvSizes[i].resize(m_mpiSurfaces[i].size());
5042 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
5043 nghbrSendSizes[i][j] = mpiSurfaceSendSize[m_mpiSurfaces[i][j]];
5044 mpiSurfaceSendSize[m_mpiSurfaces[i][j]] = 0;
5049 std::vector<MPI_Request> metaDataRecvRequests(m_noExchangeNghbrDomains);
5050 std::vector<MPI_Request> metaDataSendRequests(m_noExchangeNghbrDomains);
5051 for(
MInt nghbrIndex = 0; nghbrIndex < m_noExchangeNghbrDomains; nghbrIndex++) {
5053 m_exchangeNghbrDomains[nghbrIndex], 0, mpiComm(), &metaDataRecvRequests[nghbrIndex], AT_,
5054 "nghbrRecvSizes[nghbrIndex][0]");
5056 m_exchangeNghbrDomains[nghbrIndex], 0, mpiComm(), &metaDataSendRequests[nghbrIndex], AT_,
5057 "nghbrSendSizes[nghbrIndex][0]");
5059 MPI_Waitall(m_noExchangeNghbrDomains, &metaDataRecvRequests[0], MPI_STATUS_IGNORE, AT_);
5063 vector<vector<MFloat>> nodeVarsSendBuffers(m_noExchangeNghbrDomains);
5064 vector<vector<MFloat>> nodeVarsRecvBuffers(m_noExchangeNghbrDomains);
5065 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
5066 const MInt sendSize = accumulate(nghbrSendSizes[i].begin(), nghbrSendSizes[i].end(), 0);
5067 const MInt recvSize = accumulate(nghbrRecvSizes[i].begin(), nghbrRecvSizes[i].end(), 0);
5068 nodeVarsSendBuffers[i].resize(sendSize);
5069 nodeVarsRecvBuffers[i].resize(recvSize);
5073 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
5074 if(nghbrSendSizes[i][j] < 1) {
5077 const MInt srfcId = m_mpiSurfaces[i][j];
5081 const MFloat* data = &m_surfaces.nodeVars(srfcId, 1 - extendSide);
5082 copy(data, data + nghbrSendSizes[i][j], &nodeVarsSendBuffers[i][size]);
5083 size += nghbrSendSizes[i][j];
5085 ASSERT(size ==
static_cast<MInt>(nodeVarsSendBuffers[i].size()),
"Data size does not match buffer size.");
5087 std::vector<MPI_Request> nodeVarsRecvRequests(m_noExchangeNghbrDomains);
5088 std::vector<MPI_Request> nodeVarsSendRequests(m_noExchangeNghbrDomains);
5091 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
5092 if(nodeVarsRecvBuffers[i].size() > 0) {
5094 m_exchangeNghbrDomains[i], m_exchangeNghbrDomains[i], mpiComm(), &nodeVarsRecvRequests[i], AT_,
5095 "nodeVarsRecvBuffers[i][0]");
5097 nodeVarsRecvRequests[i] = MPI_REQUEST_NULL;
5099 if(nodeVarsSendBuffers[i].size() > 0) {
5101 m_exchangeNghbrDomains[i], domainId(), mpiComm(), &nodeVarsSendRequests[i], AT_,
5102 "nodeVarsSendBuffers[i][0]");
5104 nodeVarsSendRequests[i] = MPI_REQUEST_NULL;
5109 MPI_Waitall(m_noExchangeNghbrDomains, &nodeVarsSendRequests[0], MPI_STATUSES_IGNORE, AT_);
5113 MPI_Waitall(m_noExchangeNghbrDomains, &nodeVarsRecvRequests[0], MPI_STATUSES_IGNORE, AT_);
5115 for(
MInt nghbrIndex = 0; nghbrIndex < m_noExchangeNghbrDomains; nghbrIndex++) {
5118 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[nghbrIndex].size(); j++) {
5119 if(nghbrRecvSizes[nghbrIndex][j] < 1) {
5122 const MInt srfcId = m_mpiSurfaces[nghbrIndex][j];
5123 surfaceHasRecvd[srfcId] = 1.0;
5125 copy_n(&nodeVarsRecvBuffers[nghbrIndex][size],
5126 nghbrRecvSizes[nghbrIndex][j],
5127 &m_surfaces.nodeVars(srfcId, 1 - extendSide));
5128 copy_n(&nodeVarsRecvBuffers[nghbrIndex][size],
5129 nghbrRecvSizes[nghbrIndex][j],
5130 &m_surfaces.nodeVars(srfcId, extendSide));
5132 const MInt recvElementId = m_surfaces.nghbrElementIds(srfcId, extendSide);
5136 MBool srfcIsReverseHrefTemp =
false;
5137 for(
MInt pos = 0; pos < noSurfs; pos++) {
5138 if(hReverseSurfaceId(srfcId, pos) > -1) {
5139 srfcIsReverseHrefTemp =
true;
5142 const MBool srfcIsReverseHref = srfcIsReverseHrefTemp;
5143 MBool allSurfsReady =
true;
5144 if(srfcIsReverseHref) {
5145 const MInt mainSurfId = m_elements.surfaceIds(recvElementId, donorDir);
5146 for(
MInt pos = 0; pos < noSurfs; pos++) {
5147 if(surfaceHasRecvd[hReverseSurfaceId(mainSurfId, pos)] < 1) {
5148 allSurfsReady =
false;
5154 const MInt updateeElementId = m_surfaces.nghbrElementIds(srfcId, extendSide);
5156 if(!elementOutsideLimit[updateeElementId]) {
5157 updateNodeVariablesSingleElement(recvElementId, extendDir, hForwardSurfaceId, hReverseSurfaceId,
5158 elementUpdated, mpiSurfaceSendSize, dummy);
5159 elementUpdated[recvElementId] =
true;
5162 size += nghbrRecvSizes[nghbrIndex][j];
5164 ASSERT(size ==
static_cast<MInt>(nodeVarsRecvBuffers[nghbrIndex].size()),
5165 "Data size does not match buffer size.");
5169 MPI_Waitall(m_noExchangeNghbrDomains, &metaDataSendRequests[0], MPI_STATUS_IGNORE, AT_);
5171 m_log <<
"done" << std::endl;
5172 updateNodeVariables();
5185template <MInt nDim,
class SysEqn>
5187 const MInt extendDir,
5190 std::vector<MBool>& elementUpdated,
5191 std::vector<MInt>& mpiSurfaceSendSize,
5192 MBool& noMoreUpdates) {
5194 const MInt*
const noNodes1D = &m_elements.noNodes1D(0);
5195 const MInt*
const surfaceNoNodes1D = &m_surfaces.noNodes1D(0);
5196 const MInt extendSide = extendDir % 2;
5197 const MInt extendDimension = (extendDir - extendSide) / 2;
5198 const MInt donorDir = 2 * extendDimension + (1 - extendSide);
5199 const MInt noVars = max(SysEqn::noNodeVars(), 1);
5200 const MInt noSurfs = 2 * (nDim - 1);
5201 const MInt maxNoNodes1D = noNodes1D[elementId];
5202 const MInt maxNoNodes1D3 = (nDim == 3) ? maxNoNodes1D : 1;
5203 MInt extendDomainId;
5204 MPI_Comm_rank(mpiComm(), &extendDomainId);
5206 const MInt donorSrfcId = m_elements.surfaceIds(elementId, donorDir);
5207 const MInt recvSrfcId = m_elements.surfaceIds(elementId, extendDir);
5208 MFloatTensor nodeVars(&m_elements.nodeVars(elementId), maxNoNodes1D, maxNoNodes1D, maxNoNodes1D3, noVars);
5209 const MInt noNodesSecondSurfDim = extendDimension == 2 ? maxNoNodes1D : maxNoNodes1D3;
5210 MFloatTensor elemNodeVarsSlice(maxNoNodes1D, noNodesSecondSurfDim, noVars);
5213 MFloatTensor recvSurfaceNodeVarsMSide(&m_surfaces.nodeVars(recvSrfcId, extendSide), maxNoNodes1D,
5214 noNodesSecondSurfDim, noVars);
5215 MFloatTensor recvSurfaceNodeVarsPSide(&m_surfaces.nodeVars(recvSrfcId, 1 - extendSide), maxNoNodes1D,
5216 noNodesSecondSurfDim, noVars);
5219 const MBool pRefinedDonorSide = m_surfaces.polyDeg(donorSrfcId) > m_elements.polyDeg(elementId);
5220 const MBool pRefinedRecvSide = m_surfaces.polyDeg(recvSrfcId) > m_elements.polyDeg(elementId);
5221 const MBool srfcIsForwardHref = hForwardSurfaceId(donorSrfcId, 0) > -1;
5225 m_sysEqn.getDefaultNodeVars(defaultNodeVars.
getPointer());
5228 for(
MInt iVars = 0; iVars < noVars; iVars++) {
5229 extendNodeVar[iVars] = m_sysEqn.extendNodeVar(iVars);
5234 MBool srfcIsReverseHrefTemp =
false;
5235 for(
MInt pos = 0; pos < noSurfs; pos++) {
5236 if(hReverseSurfaceId(donorSrfcId, pos) > -1) {
5237 srfcIsReverseHrefTemp =
true;
5240 const MBool srfcIsReverseHref = srfcIsReverseHrefTemp;
5241 MBool srfcIsMainPosTemp =
false;
5242 for(
MInt pos = 0; pos < noSurfs; pos++) {
5243 if(hReverseSurfaceId(donorSrfcId, pos) == donorSrfcId) {
5244 srfcIsMainPosTemp =
true;
5247 const MBool srfcIsMainPos = srfcIsMainPosTemp;
5248 if(srfcIsForwardHref) {
5250 for(
MInt pos = 1; pos < noSurfs; pos++) {
5251 const MInt hSrfcId = hForwardSurfaceId(donorSrfcId, pos);
5252 const MInt hElemId = m_surfaces.nghbrElementIds(hSrfcId, extendSide);
5253 updateNodeVariablesSingleElement(hElemId, extendDir, hForwardSurfaceId, hReverseSurfaceId, elementUpdated,
5254 mpiSurfaceSendSize, noMoreUpdates);
5258 if(!srfcIsReverseHref) {
5260 if(pRefinedDonorSide) {
5261 const MInt maxNoNodes1DFine = surfaceNoNodes1D[donorSrfcId];
5262 const MInt maxNoNodes1D3Fine = nDim == 3 ? maxNoNodes1DFine : 1;
5263 const MInt noNodesSecondSurfDimFine = extendDimension == 2 ? maxNoNodes1DFine : maxNoNodes1D3Fine;
5266 MFloatTensor projectedRight(maxNoNodes1DFine, noNodesSecondSurfDimFine, noVars);
5267 calcMortarProjection<dg::mortar::reverse, SysEqn::noNodeVars()>(donorSrfcId, donorDir,
5268 &m_surfaces.nodeVars(donorSrfcId, 1 - extendSide),
5269 &projectedRight[0], m_elements, m_surfaces);
5270 for(
MInt a = 0;
a < maxNoNodes1D; ++
a) {
5271 for(
MInt b = 0;
b < maxNoNodes1D3; ++
b) {
5272 for(
MInt iVars = 0; iVars < noVars; iVars++) {
5273 elemNodeVarsSlice(
a,
b, iVars) = projectedRight(
a,
b, iVars);
5278 copy_n(&m_surfaces.nodeVars(donorSrfcId, 1 - extendSide),
5279 maxNoNodes1D * noNodesSecondSurfDim * noVars,
5280 &elemNodeVarsSlice[0]);
5283 for(
MInt a = 0;
a < maxNoNodes1D; ++
a) {
5284 for(
MInt b = 0;
b < maxNoNodes1D; ++
b) {
5285 for(
MInt c = 0; c < maxNoNodes1D3; ++c) {
5287 const MInt iNodes2D = extendDimension == 0 ?
b :
a;
5288 const MInt iNodes3D = extendDimension == 2 ?
b : c;
5289 for(
MInt iVars = 0; iVars < noVars; iVars++) {
5290 nodeVars(
a,
b, c, iVars) =
5291 (extendNodeVar[iVars]) ? elemNodeVarsSlice(iNodes2D, iNodes3D, iVars) : defaultNodeVars(iVars);
5297 elementUpdated[elementId] =
true;
5298 noMoreUpdates =
false;
5300 if(pRefinedRecvSide) {
5303 if(hForwardSurfaceId(recvSrfcId, 0) <= 0) {
5304 const MInt maxNoNodes1DSurf = surfaceNoNodes1D[recvSrfcId];
5305 const MInt maxNoNodes1D3Surf = nDim == 3 ? maxNoNodes1DSurf : 1;
5306 const MInt noNodesSecondRefSurfDim = extendDimension == 2 ? maxNoNodes1DSurf : maxNoNodes1D3Surf;
5307 MFloatTensor projected(maxNoNodes1DSurf, noNodesSecondRefSurfDim, noVars);
5309 calcMortarProjection<dg::mortar::forward, SysEqn::noNodeVars()>(
5310 recvSrfcId, extendDir, &m_surfaces.nodeVars(donorSrfcId, 1 - extendSide), &projected[0], m_elements,
5313 copy_n(&projected[0],
5314 maxNoNodes1DSurf * noNodesSecondRefSurfDim * noVars,
5315 &m_surfaces.nodeVars(recvSrfcId, extendSide));
5316 copy_n(&projected[0],
5317 maxNoNodes1DSurf * noNodesSecondRefSurfDim * noVars,
5318 &m_surfaces.nodeVars(recvSrfcId, 1 - extendSide));
5320 if(recvSrfcId >= m_mpiSurfacesOffset) {
5321 mpiSurfaceSendSize[recvSrfcId] = maxNoNodes1DSurf * maxNoNodes1D3Surf * SysEqn::noNodeVars();
5325 copy_n(&elemNodeVarsSlice[0], maxNoNodes1D * noNodesSecondSurfDim * noVars, &recvSurfaceNodeVarsMSide[0]);
5326 copy_n(&elemNodeVarsSlice[0], maxNoNodes1D * noNodesSecondSurfDim * noVars, &recvSurfaceNodeVarsPSide[0]);
5328 if(recvSrfcId >= m_mpiSurfacesOffset) {
5329 mpiSurfaceSendSize[recvSrfcId] = maxNoNodes1D * maxNoNodes1D3 * SysEqn::noNodeVars();
5333 if(srfcIsReverseHref && srfcIsMainPos) {
5335 MBool allDonorsHaveUpdatedMeans =
true;
5336 for(
MInt pos = 0; pos < noSurfs; pos++) {
5337 const MInt hElemId = m_surfaces.nghbrElementIds(hReverseSurfaceId(donorSrfcId, pos), 1 - extendSide);
5341 allDonorsHaveUpdatedMeans =
true;
5344 if(elementUpdated[hElemId] ==
false) {
5345 allDonorsHaveUpdatedMeans =
false;
5346 elementUpdated[elementId] =
false;
5347 noMoreUpdates =
false;
5351 if(allDonorsHaveUpdatedMeans) {
5352 const MInt maxNoNodes1DSurf = surfaceNoNodes1D[donorSrfcId];
5353 const MInt maxNoNodes1D3Surf = nDim == 3 ? maxNoNodes1DSurf : 1;
5354 const MInt noNodesSecondRefSurfDim = extendDimension == 2 ? maxNoNodes1DSurf : maxNoNodes1D3Surf;
5355 const MInt coarseNextSrfcId = m_elements.surfaceIds(elementId, extendDir);
5357 elementUpdated[elementId] =
true;
5359 if(coarseNextSrfcId >= m_mpiSurfacesOffset) {
5360 mpiSurfaceSendSize[coarseNextSrfcId] = maxNoNodes1D * maxNoNodes1D3 * SysEqn::noNodeVars();
5362 noMoreUpdates =
false;
5363 MFloatTensor projectedSum(maxNoNodes1DSurf, noNodesSecondRefSurfDim, noVars);
5364 projectedSum.
set(0.0);
5365 for(
MInt pos = 0; pos < noSurfs; pos++) {
5366 const MInt hSrfId = hReverseSurfaceId(donorSrfcId, pos);
5367 MFloatTensor projected(maxNoNodes1DSurf, noNodesSecondRefSurfDim, noVars);
5369 calcMortarProjection<dg::mortar::reverse, SysEqn::noNodeVars()>(
5370 hSrfId, donorDir, &m_surfaces.nodeVars(hSrfId, 1 - extendSide), &projected[0], m_elements, m_surfaces);
5371 for(
MInt i = 0; i < maxNoNodes1D; i++) {
5372 for(
MInt j = 0; j < noNodesSecondSurfDim; j++) {
5373 for(
MInt k = 0; k < noVars; k++) {
5374 projectedSum(i, j, k) += projected(i, j, k);
5383 MFloatTensor hNodeVarsLeft(&m_elements.nodeVars(elementId), maxNoNodes1D, maxNoNodes1D, maxNoNodes1D3, noVars);
5384 MFloatTensor hNodeVarsLeftSurf1(&m_surfaces.nodeVars(coarseNextSrfcId, extendSide), maxNoNodes1D,
5385 noNodesSecondSurfDim, noVars);
5386 MFloatTensor hNodeVarsLeftSurf2(&m_surfaces.nodeVars(coarseNextSrfcId, 1 - extendSide), maxNoNodes1D,
5387 noNodesSecondSurfDim, noVars);
5388 for(
MInt a = 0;
a < maxNoNodes1D; ++
a) {
5389 for(
MInt b = 0;
b < maxNoNodes1D; ++
b) {
5390 for(
MInt c = 0; c < maxNoNodes1D3; ++c) {
5391 const MInt iNodes2D = extendDimension == 0 ?
b :
a;
5392 const MInt iNodes3D = extendDimension == 2 ?
b : c;
5393 for(
MInt iVars = 0; iVars < noVars; iVars++) {
5394 if(extendNodeVar[iVars]) {
5395 hNodeVarsLeft(
a,
b, c, iVars) = projectedSum(iNodes2D, iNodes3D, iVars);
5396 elemNodeVarsSlice(iNodes2D, iNodes3D, iVars) = projectedSum(iNodes2D, iNodes3D, iVars);
5397 hNodeVarsLeftSurf1(iNodes2D, iNodes3D, iVars) = projectedSum(iNodes2D, iNodes3D, iVars);
5398 hNodeVarsLeftSurf2(iNodes2D, iNodes3D, iVars) = projectedSum(iNodes2D, iNodes3D, iVars);
5400 hNodeVarsLeft(
a,
b, c, iVars) = defaultNodeVars(iVars);
5401 elemNodeVarsSlice(iNodes2D, iNodes3D, iVars) = defaultNodeVars(iVars);
5402 hNodeVarsLeftSurf1(iNodes2D, iNodes3D, iVars) = defaultNodeVars(iVars);
5403 hNodeVarsLeftSurf2(iNodes2D, iNodes3D, iVars) = defaultNodeVars(iVars);
5409 elementUpdated[elementId] =
true;
5410 noMoreUpdates =
false;
5416 if(hForwardSurfaceId(recvSrfcId, 0) > 0) {
5417 const MInt maxNoNodes1DFine = surfaceNoNodes1D[recvSrfcId];
5418 const MInt maxNoNodes1D3Fine = nDim == 3 ? maxNoNodes1DFine : 1;
5419 const MInt noNodesSecondFineSurfDim = extendDimension == 2 ? maxNoNodes1DFine : maxNoNodes1D3Fine;
5420 for(
MInt pos = 0; pos < noSurfs; pos++) {
5421 const MInt hSrfId = hForwardSurfaceId(recvSrfcId, pos);
5422 copy_n(&elemNodeVarsSlice(0, 0, 0),
5423 maxNoNodes1D * noNodesSecondSurfDim * noVars,
5424 &m_surfaces.nodeVars(hSrfId, 1 - extendSide));
5426 for(
MInt pos = 0; pos < noSurfs; pos++) {
5427 const MInt hSrfId = hForwardSurfaceId(recvSrfcId, pos);
5428 MFloatTensor projected(maxNoNodes1DFine, noNodesSecondFineSurfDim, noVars);
5430 calcMortarProjection<dg::mortar::forward, SysEqn::noNodeVars()>(
5431 hSrfId, extendDir, &m_surfaces.nodeVars(hSrfId, 1 - extendSide), &projected[0], m_elements, m_surfaces);
5432 MFloatTensor nodeVarsNextSrfPSide(&m_surfaces.nodeVars(hSrfId, extendSide), maxNoNodes1DFine,
5433 noNodesSecondFineSurfDim, noVars);
5434 MFloatTensor nodeVarsNextSrfMSide(&m_surfaces.nodeVars(hSrfId, 1 - extendSide), maxNoNodes1DFine,
5435 noNodesSecondFineSurfDim, noVars);
5436 for(
MInt iNodes2D = 0; iNodes2D < maxNoNodes1DFine; iNodes2D++) {
5437 for(
MInt iNodes3D = 0; iNodes3D < noNodesSecondFineSurfDim; iNodes3D++) {
5438 for(
MInt iVars = 0; iVars < noVars; iVars++) {
5439 nodeVarsNextSrfPSide(iNodes2D, iNodes3D, iVars) = projected(iNodes2D, iNodes3D, iVars);
5440 nodeVarsNextSrfMSide(iNodes2D, iNodes3D, iVars) = projected(iNodes2D, iNodes3D, iVars);
5445 if(recvSrfcId >= m_mpiSurfacesOffset) {
5446 mpiSurfaceSendSize[hSrfId] = maxNoNodes1DFine * maxNoNodes1D3Fine * noVars;
5459template <MInt nDim,
class SysEqn>
5464 const MString polynomialType =
"Legendre";
5465 const MString integrationMethod = (m_dgIntegrationMethod == 0) ?
"Gauss" :
"Gauss-Lobatto";
5470 Frame summary(
" SOLVER " + std::to_string(solverId()) +
" INITIALIZATION SUMMARY AT TIME STEP "
5471 + std::to_string(m_timeStep));
5475 problem.
addData(
"System of equations", m_sysEqn.sysEqnName());
5476 problem.
addData(
"Number of dimensions", nDim);
5477 Data& vars = problem.
addData(
"Number of variables", m_sysEqn.noVars());
5478 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
5480 vars.
addData(
"Conservative variable name(s)", m_sysEqn.consVarNames(i));
5482 vars.
addData(
"", m_sysEqn.consVarNames(i));
5485 problem.
addData(
"Restart", m_restart);
5487 problem.
addData(
"Initial condition", getRestartFileName(m_timeStep, m_useNonSpecifiedRestartFile));
5489 problem.
addData(
"Initial condition", m_sysEqn.m_initialCondition);
5491 problem.
addData(
"Start time (non-dimensionalized)", m_time);
5492 problem.
addData(
"Final time (non-dimensionalized)", m_finalTime);
5493 problem.
addData(
"CFL", m_sysEqn.cfl());
5494 problem.
addData(
"Recalculation interval for time step", m_calcTimeStepInterval);
5495 problem.
addData(
"Time step", m_timeStep);
5496 problem.
addData(
"Maximum number of time steps", m_timeSteps);
5499 Group& discret = summary.
addGroup(
"DISCRETIZATION SUMMARY");
5500 discret.
addData(
"Initial polynomial degree", m_initPolyDeg);
5501 discret.
addData(
"Minimum polynomial degree (limit)", m_minPolyDeg);
5502 discret.
addData(
"Maximum polynomial degree (limit)", m_maxPolyDeg);
5503 discret.
addData(
"Polynomial type", polynomialType);
5504 discret.
addData(
"Integration method", integrationMethod);
5505 MString timeIntegrationScheme;
5506 switch(m_dgTimeIntegrationScheme) {
5509 timeIntegrationScheme =
"Carpenter 4/5";
5513 timeIntegrationScheme =
"Toulorge C 4/8";
5517 timeIntegrationScheme =
"Niegemann 4/14";
5521 timeIntegrationScheme =
"Niegemann 4/13";
5525 timeIntegrationScheme =
"Toulorge C 3/7";
5529 timeIntegrationScheme =
"Toulorge F 4/8";
5533 discret.
addData(
"Time integration scheme", timeIntegrationScheme);
5536 Group& parallel = summary.
addGroup(
"PARALLELIZATION SUMMARY");
5537 parallel.
addData(
"Domain id", domainId());
5538 parallel.
addData(
"Number of neighbor domains", grid().noNeighborDomains());
5539 parallel.
addData(
"Number of exchange neighbor domains", m_noExchangeNghbrDomains);
5540 parallel.
addData(
"Total number of domains", noDomains());
5543 Group& gridLocal = summary.
addGroup(
"GRID SUMMARY (LOCAL)");
5544 gridLocal.
addData(
"Minimum used polynomial degree", m_statLocalMinPolyDeg);
5545 gridLocal.
addData(
"Maximum used polynomial degree", m_statLocalMaxPolyDeg);
5546 gridLocal.
addData(
"Average used polynomial degree", 1.0 * m_statLocalPolyDegSum / m_statLocalNoActiveCells);
5548 gridLocal.
addData(
"Minimum grid level", m_statLocalMinLevel);
5549 gridLocal.
addData(
"Maximum grid level", m_statLocalMaxLevel);
5551 Data& noCellsLocal = gridLocal.
addData(
"Number of cells", m_statLocalNoCells);
5552 noCellsLocal.
addData(
"internal cells", m_statLocalNoInternalCells);
5553 noCellsLocal.
addData(
"halo cells", m_statLocalNoHaloCells);
5554 gridLocal.
addData(
"Number of active cells", m_statLocalNoActiveCells);
5555 gridLocal.
addData(
"Number of active DOFs", m_statLocalNoActiveDOFs);
5557 const MInt noPolyDegs = m_maxPolyDeg - m_minPolyDeg + 1;
5558 if(noPolyDegs > 1) {
5559 for(
MInt i = 0; i < noPolyDegs; i++) {
5560 MString name =
"Number of active DOFs polyDeg=" + std::to_string(m_minPolyDeg + i);
5561 gridLocal.
addData(name, m_statLocalNoActiveDOFsPolyDeg[i]);
5564 gridLocal.
addData(
"Max. number of cells (collector size)", m_statLocalMaxNoCells);
5565 gridLocal.
addData(
"Memory utilization", 100.0 * m_statLocalNoCells / m_statLocalMaxNoCells);
5567 gridLocal.
addData(
"Number of elements", m_statLocalNoElements);
5568 gridLocal.
addData(
"Number of helements", m_statLocalNoHElements);
5570 Data& surfaces = gridLocal.
addData(
"Number of surfaces", m_statLocalNoSurfaces);
5571 surfaces.
addData(
"boundary surfaces", m_statLocalNoBoundarySurfaces);
5572 surfaces.
addData(
"inner surfaces", m_statLocalNoInnerSurfaces);
5573 surfaces.
addData(
"MPI surfaces", m_statLocalNoMpiSurfaces);
5574 gridLocal.
addData(
"Max. number of surfaces (collector size)", m_statLocalMaxNoSurfaces);
5575 gridLocal.
addData(
"Memory utilization", 100.0 * m_statLocalNoSurfaces / m_statLocalMaxNoSurfaces);
5578 Group& gridGlobal = summary.
addGroup(
"GRID SUMMARY (GLOBAL)");
5579 gridGlobal.
addData(
"Minimum used polynomial degree", m_statGlobalMinPolyDeg);
5580 gridGlobal.
addData(
"Maximum used polynomial degree", m_statGlobalMaxPolyDeg);
5581 gridGlobal.
addData(
"Average used polynomial degree", 1.0 * m_statGlobalPolyDegSum / m_statGlobalNoActiveCells);
5583 gridGlobal.
addData(
"Minimum grid level", m_statGlobalMinLevel);
5584 gridGlobal.
addData(
"Maximum grid level", m_statGlobalMaxLevel);
5586 Data& noCellsGlbl = gridGlobal.
addData(
"Number of cells", m_statGlobalNoCells);
5587 noCellsGlbl.
addData(
"internal cells", m_statGlobalNoInternalCells);
5588 noCellsGlbl.
addData(
"halo cells", m_statGlobalNoHaloCells);
5589 gridGlobal.
addData(
"Number of active cells", m_statGlobalNoActiveCells);
5590 gridGlobal.
addData(
"Number of active DOFs", m_statGlobalNoActiveDOFs);
5591 gridGlobal.
addData(
"Max. number of cells (collector size)", m_statGlobalMaxNoCells);
5592 gridGlobal.
addData(
"Memory utilization", 100.0 * m_statGlobalNoCells / m_statGlobalMaxNoCells);
5594 gridGlobal.
addData(
"Number of elements", m_statGlobalNoElements);
5595 gridGlobal.
addData(
"Number of helements", m_statGlobalNoHElements);
5597 Data& surfacesGlobal = gridGlobal.
addData(
"Number of surfaces", m_statGlobalNoSurfaces);
5598 surfacesGlobal.
addData(
"boundary surfaces", m_statGlobalNoBoundarySurfaces);
5599 surfacesGlobal.
addData(
"inner surfaces", m_statGlobalNoInnerSurfaces);
5600 surfacesGlobal.
addData(
"MPI surfaces", m_statGlobalNoMpiSurfaces);
5601 surfacesGlobal.
addData(
"Surfaces marked for p-ref", m_statGlobalNoPrefSurfs);
5602 surfacesGlobal.
addData(
"Surfaces marked for h-ref", m_statGlobalNoHrefSurfs);
5603 gridGlobal.
addData(
"Max. number of surfaces (collector size)", m_statGlobalMaxNoSurfaces);
5604 gridGlobal.
addData(
"Memory utilization", 100.0 * m_statGlobalNoSurfaces / m_statGlobalMaxNoSurfaces);
5608 for(
auto&& bc : m_boundaryConditions) {
5609 boundary.
addData(bc->name(), bc->id());
5613 Group& cutOffBoundary = summary.
addGroup(
"CUT-OFF BOUNDARY CONDITIONS");
5614 cutOffBoundary.
addData(
"Enabled?", m_useCutOffBoundaries);
5615 if(m_useCutOffBoundaries) {
5616 std::array<MString, 6> dirs = {{
"-x",
"+x",
"-y",
"+y",
"-z",
"+z"}};
5617 for(
MInt i = 0; i < 2 * nDim; i++) {
5618 const MString bc = dirs[i] +
" (bcId: " + to_string(m_cutOffBoundaryConditionIds[i]) +
")";
5619 cutOffBoundary.
addData(bc, m_statGlobalNoCutOffBoundarySurfaces[i]);
5625 sbp.
addData(
"Enabled?", m_sbpMode);
5626 sbp.
addData(
"Operator", m_sbpOperator);
5633 cout << s << std::endl;
5648template <MInt nDim,
class SysEqn>
5653 ss << setw(8) << setfill(
'0') << m_timeStep;
5654 saveSolutionFile(ss.str());
5667template <MInt nDim,
class SysEqn>
5670 RECORD_TIMER_START(m_timers[Timers::Accumulated]);
5671 RECORD_TIMER_START(m_timers[Timers::IO]);
5672 RECORD_TIMER_START(m_timers[Timers::SaveSolutionFile]);
5674 m_log <<
"Saving solution file ... ";
5677 using namespace parallel_io;
5679 outputDir() +
"solution_" + getIdentifier(
g_multiSolverGrid,
"b") + suffix + ParallelIo::fileExt();
5680 const MInt noElements = m_elements.size();
5681 ParallelIo parallelIo(fileName, PIO_REPLACE, mpiComm());
5684 if(hasAdaptivePref()) {
5685 parallelIo.setAttribute(suffix + grid().gridInputFileName(),
"gridFile");
5687 parallelIo.setAttribute(grid().gridInputFileName(),
"gridFile");
5691 parallelIo.setAttribute(solverId(),
"solverId");
5696 for(
MInt cellId = 0, offset = 0, elementId = 0; elementId < noElements; cellId++) {
5699 if(cellId != m_elements.cellId(elementId)) {
5700 offset +=
ipow(m_initNoNodes1D, nDim);
5705 elementOffset[elementId] = offset;
5708 const MInt noNodesXD = m_elements.noNodesXD(elementId);
5709 offset += noNodesXD;
5712 elementNodes[elementId] = noNodesXD;
5719 const MInt localNoNodes = m_noTotalNodesXD + (grid().noInternalCells() - noElements) *
ipow(m_initNoNodes1D, nDim);
5722 ParallelIo::size_type nodesOffset, globalNoNodes;
5723 ParallelIo::calcOffset(localNoNodes, &nodesOffset, &globalNoNodes, mpiComm());
5726 parallelIo.setAttribute(
"DG",
"solverType");
5728 parallelIo.setAttribute((
MInt)m_sbpMode,
"sbpMode");
5730 parallelIo.setAttribute(m_timeStep,
"timeStep");
5731 parallelIo.setAttribute(m_time,
"time");
5735 parallelIo.defineArray(PIO_UCHAR,
"polyDegs", grid().domainOffset(noDomains()));
5738 parallelIo.defineArray(PIO_UCHAR,
"noNodes1D", grid().domainOffset(noDomains()));
5743 MString dgIntegrationMethod = m_sbpMode ?
"DG_INTEGRATE_GAUSS_LOBATTO" :
"DG_INTEGRATE_GAUSS";
5744 dgIntegrationMethod =
5745 Context::getSolverProperty<MString>(
"dgIntegrationMethod", m_solverId, AT_, &dgIntegrationMethod);
5746 MString dgPolynomialType =
"DG_POLY_LEGENDRE";
5747 dgPolynomialType = Context::getSolverProperty<MString>(
"dgPolynomialType", m_solverId, AT_, &dgPolynomialType);
5750 parallelIo.setAttribute(dgIntegrationMethod,
"dgIntegrationMethod",
"polyDegs");
5751 parallelIo.setAttribute(dgPolynomialType,
"dgPolynomialType",
"polyDegs");
5757 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
5758 const MString name =
"variables" + to_string(varId++);
5759 parallelIo.defineArray(PIO_FLOAT, name, globalNoNodes);
5760 parallelIo.setAttribute(m_sysEqn.consVarNames(i),
"name", name);
5764 if(m_saveNodeVariablesToSolutionFile) {
5765 for(
MInt i = 0; i < SysEqn::noNodeVars(); i++) {
5766 const MString name =
"variables" + to_string(varId++);
5767 parallelIo.defineArray(PIO_FLOAT, name, globalNoNodes);
5768 parallelIo.setAttribute(m_sysEqn.nodeVarNames(i),
"name", name);
5773 if(m_writeTimeDerivative) {
5774 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
5775 const MString name =
"variables" + to_string(varId++);
5776 parallelIo.defineArray(PIO_FLOAT, name, globalNoNodes);
5777 parallelIo.setAttribute(m_sysEqn.consVarNames(i) +
"_tDeriv",
"name", name);
5781 if(useSponge() && m_writeSpongeEta > 0) {
5782 const MString name =
"variables" + to_string(varId);
5783 parallelIo.defineArray(PIO_FLOAT, name, globalNoNodes);
5784 parallelIo.setAttribute(
"spongeEta",
"name", name);
5789 parallelIo.setOffset(grid().noInternalCells(), grid().domainOffset(domainId()));
5791 fill_n(&polyDegs[0], grid().noInternalCells(), m_initPolyDeg);
5792 for(
MInt elementId = 0; elementId < noElements; elementId++) {
5793 const MInt cellId = m_elements.cellId(elementId);
5794 polyDegs.
p[cellId] = m_elements.polyDeg(elementId);
5796 parallelIo.writeArray(polyDegs.
begin(),
"polyDegs");
5800 parallelIo.setOffset(grid().noInternalCells(), grid().domainOffset(domainId()));
5802 fill_n(&noNodes1D[0], grid().noInternalCells(), m_initNoNodes1D);
5803 for(
MInt elementId = 0; elementId < noElements; elementId++) {
5804 const MInt cellId = m_elements.cellId(elementId);
5805 noNodes1D.
p[cellId] = m_elements.noNodes1D(elementId);
5807 parallelIo.writeArray(noNodes1D.
begin(),
"noNodes1D");
5811 parallelIo.setOffset(localNoNodes, nodesOffset);
5814 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
5816 fill(buffer.
begin(), buffer.
end(), 0.0);
5817 for(
MInt e = 0; e < noElements; e++) {
5818 MFloat*
const b = &buffer[elementOffset[e]];
5819 const MFloat*
const v = &m_elements.variables(e) + i;
5820 for(
MInt n = 0; n < elementNodes[e]; n++) {
5821 b[n] = v[n * m_sysEqn.noVars()];
5824 const MString name =
"variables" + to_string(varId++);
5825 parallelIo.writeArray(&buffer[0], name);
5829 if(m_saveNodeVariablesToSolutionFile) {
5830 for(
MInt i = 0; i < SysEqn::noNodeVars(); i++) {
5832 fill(buffer.
begin(), buffer.
end(), 0.0);
5833 for(
MInt e = 0; e < noElements; e++) {
5834 MFloat*
const b = &buffer[elementOffset[e]];
5835 const MFloat*
const v = &m_elements.nodeVars(e) + i;
5836 for(
MInt n = 0; n < elementNodes[e]; n++) {
5837 b[n] = v[n * SysEqn::noNodeVars()];
5840 const MString name =
"variables" + to_string(varId++);
5841 parallelIo.writeArray(&buffer[0], name);
5846 if(m_writeTimeDerivative) {
5847 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
5849 fill(buffer.
begin(), buffer.
end(), 0.0);
5850 for(
MInt e = 0; e < noElements; e++) {
5851 MFloat*
const b = &buffer[elementOffset[e]];
5852 const MFloat*
const r = &m_elements.rightHandSide(e) + i;
5853 for(
MInt n = 0; n < elementNodes[e]; n++) {
5854 b[n] = r[n * m_sysEqn.noVars()];
5857 const MString name =
"variables" + to_string(varId++);
5858 parallelIo.writeArray(&buffer[0], name);
5863 if(useSponge() && m_writeSpongeEta > 0) {
5864 const MInt noSpongeElements = sponge().noSpongeElements();
5866 fill(buffer.
begin(), buffer.
end(), 0.0);
5867 for(
MInt e = 0; e < noSpongeElements; e++) {
5868 MInt elementId = sponge().elementId(e);
5869 MFloat*
const b = &buffer[elementOffset[elementId]];
5870 for(
MInt n = 0; n < elementNodes[elementId]; n++) {
5871 b[n] = sponge().spongeEta(e, n);
5874 const MString name =
"variables" + to_string(varId);
5875 parallelIo.writeArray(&buffer[0], name);
5877 m_log <<
"done" << endl;
5879 RECORD_TIMER_STOP(m_timers[Timers::SaveSolutionFile]);
5880 RECORD_TIMER_STOP(m_timers[Timers::IO]);
5881 RECORD_TIMER_STOP(m_timers[Timers::Accumulated]);
5892template <MInt nDim,
class SysEqn>
5895 RECORD_TIMER_START(m_timers[Timers::Accumulated]);
5896 RECORD_TIMER_START(m_timers[Timers::IO]);
5897 RECORD_TIMER_START(m_timers[Timers::SaveRestartFile]);
5900 const MString fileName = getRestartFileName(m_timeStep, m_useNonSpecifiedRestartFile);
5903 ss <<
"restart_" << getIdentifier(
g_multiSolverGrid,
"b") << setw(8) << setfill(
'0') << m_timeStep;
5904 const MString suffix = ss.str();
5907 const MInt noElements = m_elements.size();
5910 for(
MInt cellId = 0, offset = 0, elementId = 0; elementId < noElements; cellId++) {
5913 if(cellId != m_elements.cellId(elementId)) {
5914 offset +=
ipow(m_initNoNodes1D, nDim);
5919 elementOffset[elementId] = offset;
5922 const MInt noNodesXD = m_elements.noNodesXD(elementId);
5923 offset += noNodesXD;
5926 elementNodes[elementId] = noNodesXD;
5933 const MInt localNoNodes = m_noTotalNodesXD + (grid().noInternalCells() - noElements) *
ipow(m_initNoNodes1D, nDim);
5936 using namespace parallel_io;
5937 ParallelIo parallelIo(fileName, PIO_REPLACE, mpiComm());
5941 if(hasAdaptivePref()) {
5942 parallelIo.setAttribute(suffix + grid().gridInputFileName(),
"gridFile");
5944 parallelIo.setAttribute(grid().gridInputFileName(),
"gridFile");
5946 parallelIo.setAttribute(
"DG",
"solverType");
5948 parallelIo.setAttribute((
MInt)m_sbpMode,
"sbpMode");
5950 parallelIo.setAttribute(solverId(),
"solverId");
5953 parallelIo.defineScalar(PIO_INT,
"timeStep");
5954 parallelIo.defineScalar(PIO_FLOAT,
"time");
5955 parallelIo.defineScalar(PIO_FLOAT,
"dt");
5958 ParallelIo::size_type nodesOffset, globalNoNodes;
5959 ParallelIo::calcOffset(localNoNodes, &nodesOffset, &globalNoNodes, mpiComm());
5963 parallelIo.defineArray(PIO_UCHAR,
"polyDegs", grid().domainOffset(noDomains()));
5966 parallelIo.defineArray(PIO_UCHAR,
"noNodes1D", grid().domainOffset(noDomains()));
5970 MString dgIntegrationMethod = m_sbpMode ?
"DG_INTEGRATE_GAUSS_LOBATTO" :
"DG_INTEGRATE_GAUSS";
5971 dgIntegrationMethod =
5972 Context::getSolverProperty<MString>(
"dgIntegrationMethod", m_solverId, AT_, &dgIntegrationMethod);
5973 MString dgPolynomialType =
"DG_POLY_LEGENDRE";
5974 dgPolynomialType = Context::getSolverProperty<MString>(
"dgPolynomialType", m_solverId, AT_, &dgPolynomialType);
5976 parallelIo.setAttribute(dgIntegrationMethod,
"dgIntegrationMethod",
"polyDegs");
5977 parallelIo.setAttribute(dgPolynomialType,
"dgPolynomialType",
"polyDegs");
5983 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
5984 const MString name =
"variables" + to_string(varId++);
5985 parallelIo.defineArray(PIO_FLOAT, name, globalNoNodes);
5986 parallelIo.setAttribute(m_sysEqn.consVarNames(i),
"name", name);
5990 if(SysEqn::hasTimeDependentNodeVars()) {
5991 for(
MInt i = 0; i < SysEqn::noNodeVars(); i++) {
5992 const MString name =
"variables" + to_string(varId++);
5993 parallelIo.defineArray(PIO_FLOAT, name, globalNoNodes);
5994 parallelIo.setAttribute(m_sysEqn.nodeVarNames(i),
"name", name);
6016 const MInt noBcIds = m_boundaryConditions.size();
6017 std::vector<ParallelIo::size_type> bcNodesOffset(noBcIds, 0), bcLocalNoNodes(noBcIds, 0), bcGlobalNoNodes(noBcIds, 0);
6019 for(
auto&& bc : m_boundaryConditions) {
6020 const MInt bcNoVars = bc->noRestartVars();
6022 bcLocalNoNodes[bcId] = bc->getLocalNoNodes();
6023 ParallelIo::calcOffset(bcLocalNoNodes[bcId], &bcNodesOffset[bcId], &bcGlobalNoNodes[bcId], mpiComm());
6024 for(
MInt i = 0; i < bcNoVars; i++) {
6025 const MString name =
"variables" + to_string(varId++);
6026 parallelIo.defineArray(PIO_FLOAT, name, bcGlobalNoNodes[bcId]);
6027 parallelIo.setAttribute(bc->restartVarName(i),
"name", name);
6036 parallelIo.writeScalar(m_timeStep,
"timeStep");
6037 parallelIo.writeScalar(m_time,
"time");
6038 parallelIo.writeScalar(m_dt,
"dt");
6041 parallelIo.setOffset(grid().noInternalCells(), grid().domainOffset(domainId()));
6043 fill_n(&polyDegs[0], grid().noInternalCells(), m_initPolyDeg);
6044 for(
MInt elementId = 0; elementId < noElements; elementId++) {
6045 const MInt cellId = m_elements.cellId(elementId);
6046 polyDegs.
p[cellId] = m_elements.polyDeg(elementId);
6048 parallelIo.writeArray(polyDegs.
begin(),
"polyDegs");
6052 parallelIo.setOffset(grid().noInternalCells(), grid().domainOffset(domainId()));
6054 fill_n(&noNodes1D[0], grid().noInternalCells(), m_initNoNodes1D);
6055 for(
MInt elementId = 0; elementId < noElements; elementId++) {
6056 const MInt cellId = m_elements.cellId(elementId);
6057 noNodes1D.
p[cellId] = m_elements.noNodes1D(elementId);
6059 parallelIo.writeArray(noNodes1D.
begin(),
"noNodes1D");
6063 parallelIo.setOffset(localNoNodes, nodesOffset);
6066 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
6069 fill(buffer.
begin(), buffer.
end(), 0.0);
6070 for(
MInt e = 0; e < noElements; e++) {
6071 MFloat*
const b = &buffer[elementOffset[e]];
6072 const MFloat*
const v = &m_elements.variables(e) + i;
6073 for(
MInt n = 0; n < elementNodes[e]; n++) {
6074 b[n] = v[n * m_sysEqn.noVars()];
6077 const MString name =
"variables" + to_string(varId++);
6078 parallelIo.writeArray(&buffer[0], name);
6082 if(SysEqn::hasTimeDependentNodeVars()) {
6083 for(
MInt i = 0; i < SysEqn::noNodeVars(); i++) {
6085 fill(buffer.
begin(), buffer.
end(), 0.0);
6086 for(
MInt e = 0; e < noElements; e++) {
6087 MFloat*
const b = &buffer[elementOffset[e]];
6088 const MFloat*
const v = &m_elements.nodeVars(e) + i;
6089 for(
MInt n = 0; n < elementNodes[e]; n++) {
6090 b[n] = v[n * SysEqn::noNodeVars()];
6093 const MString name =
"variables" + to_string(varId++);
6094 parallelIo.writeArray(&buffer[0], name);
6144 for(
auto&& bc : m_boundaryConditions) {
6145 parallelIo.setOffset(bcLocalNoNodes[bcId], bcNodesOffset[bcId]);
6146 for(
MInt i = 0; i < bc->noRestartVars(); i++) {
6149 fill(buffer.
begin(), buffer.
end(), 0.0);
6152 bc->getRestartVariable(i, &buffer[0]);
6153 const MString name =
"variables" + to_string(varId++);
6154 parallelIo.writeArray(&buffer[0], name);
6159 RECORD_TIMER_STOP(m_timers[Timers::SaveRestartFile]);
6160 RECORD_TIMER_STOP(m_timers[Timers::IO]);
6161 RECORD_TIMER_STOP(m_timers[Timers::Accumulated]);
6172template <MInt nDim,
class SysEqn>
6176 m_log <<
"Loading restart file... ";
6179 const MString fileName = getRestartFileName(m_restartTimeStep, m_useNonSpecifiedRestartFile);
6182 using namespace parallel_io;
6183 ParallelIo parallelIo(fileName, PIO_READ, mpiComm());
6187 parallelIo.getAttribute(&gridFileName,
"gridFile");
6190 parallelIo.readScalar(&m_timeStep,
"timeStep");
6191 parallelIo.readScalar(&m_time,
"time");
6192 parallelIo.readScalar(&m_dt,
"dt");
6196 if(!parallelIo.hasDataset(
"polyDegs", 1)) {
6197 TERMM(1,
"ERROR: restart file is not valid (missing polynomial degree data).");
6199 parallelIo.setOffset(grid().noInternalCells(), grid().domainOffset(domainId()));
6200 parallelIo.readArray(&polyDegs[0],
"polyDegs");
6203 const MInt noElements = m_elements.size();
6204 m_noTotalNodesXD = 0;
6205 for(
MInt elementId = 0; elementId < noElements; elementId++) {
6206 const MInt cellId = m_elements.cellId(elementId);
6207 const MInt polyDeg = polyDegs[cellId];
6208 const MInt noNodes1D = polyDeg + 1;
6209 const MInt noNodesXD =
ipow(polyDeg + 1, nDim);
6210 m_elements.polyDeg(elementId) = polyDeg;
6211 m_elements.noNodes1D(elementId) = noNodes1D;
6212 m_noTotalNodesXD += noNodesXD;
6217 for(
MInt cellId = 0, offset = 0, elementId = 0; elementId < noElements; cellId++) {
6220 if(cellId != m_elements.cellId(elementId)) {
6221 offset +=
ipow(m_initNoNodes1D, nDim);
6226 elementOffset[elementId] = offset;
6229 const MInt noNodesXD = m_elements.noNodesXD(elementId);
6230 offset += noNodesXD;
6240 const MInt localNoNodes = m_noTotalNodesXD + (grid().noInternalCells() - noElements) *
ipow(m_initNoNodes1D, nDim);
6243 ParallelIo::size_type nodesOffset, globalNoNodes;
6244 ParallelIo::calcOffset(localNoNodes, &nodesOffset, &globalNoNodes, mpiComm());
6247 parallelIo.setOffset(localNoNodes, nodesOffset);
6251 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
6253 const MString name =
"variables" + to_string(varId++);
6254 parallelIo.readArray(&buffer[0], name);
6255 for(
MInt e = 0; e < noElements; e++) {
6256 const MFloat*
const b = &buffer[elementOffset[e]];
6257 MFloat*
const v = &m_elements.variables(e) + i;
6258 const MInt noNodesXD = m_elements.noNodesXD(e);
6259 for(
MInt n = 0; n < noNodesXD; n++) {
6260 v[n * m_sysEqn.noVars()] =
b[n];
6266 if(SysEqn::hasTimeDependentNodeVars()) {
6267 for(
MInt i = 0; i < m_sysEqn.noNodeVars(); i++) {
6269 const MString name =
"variables" + to_string(varId++);
6270 parallelIo.readArray(&buffer[0], name);
6271 for(
MInt e = 0; e < noElements; e++) {
6272 const MFloat*
const b = &buffer[elementOffset[e]];
6273 MFloat*
const v = &m_elements.nodeVars(e) + i;
6274 const MInt noNodesXD = m_elements.noNodesXD(e);
6275 for(
MInt n = 0; n < noNodesXD; n++) {
6276 v[n * m_sysEqn.noNodeVars()] =
b[n];
6320 for(
auto&& bc : m_boundaryConditions) {
6321 const MInt bcNoVars = bc->noRestartVars();
6323 const MInt bcLocalNoNodes = bc->getLocalNoNodes();
6324 ParallelIo::size_type bcNodesOffset, bcGlobalNoNodes;
6325 ParallelIo::calcOffset(bcLocalNoNodes, &bcNodesOffset, &bcGlobalNoNodes, mpiComm());
6326 parallelIo.setOffset(bcLocalNoNodes, bcNodesOffset);
6327 for(
MInt i = 0; i < bcNoVars; i++) {
6330 const MString name =
"variables" + to_string(varId++);
6331 parallelIo.readArray(&buffer[0], name);
6334 bc->setRestartVariable(i, &buffer[0]);
6339 m_log <<
"done" << endl;
6360template <MInt nDim,
class SysEqn>
6363 const vector<MString>& varNames,
6364 const MFloat*
const data)
const {
6370 TERMM(1,
"noVars must by >= 1");
6374 const MInt noElements = m_elements.size();
6377 for(
MInt cellId = 0, offset = 0, elementId = 0; elementId < noElements; cellId++) {
6380 if(cellId != m_elements.cellId(elementId)) {
6381 offset +=
ipow(m_initNoNodes1D, nDim);
6386 elementOffset[elementId] = offset;
6389 const MInt noNodesXD = m_elements.noNodesXD(elementId);
6390 offset += noNodesXD;
6393 elementNodes[elementId] = noNodesXD;
6400 const MInt localNoNodes = m_noTotalNodesXD + (grid().noInternalCells() - noElements) *
ipow(m_initNoNodes1D, nDim);
6403 const MString fileName = outputDir() + fileNameBase + ParallelIo::fileExt();
6404 using namespace parallel_io;
6405 ParallelIo parallelIo(fileName, PIO_REPLACE, mpiComm());
6408 parallelIo.setAttribute(grid().gridInputFileName(),
"gridFile");
6409 parallelIo.setAttribute(
"DG",
"solverType");
6411 parallelIo.setAttribute((
MInt)m_sbpMode,
"sbpMode");
6415 if(grid().raw().treeb().
noSolvers() > 1) {
6416 parallelIo.setAttribute(solverId(),
"solverId");
6420 ParallelIo::size_type nodesOffset, globalNoNodes;
6421 ParallelIo::calcOffset(localNoNodes, &nodesOffset, &globalNoNodes, mpiComm());
6425 parallelIo.defineArray(PIO_UCHAR,
"polyDegs", grid().domainOffset(noDomains()));
6428 parallelIo.defineArray(PIO_UCHAR,
"noNodes1D", grid().domainOffset(noDomains()));
6432 MString dgIntegrationMethod = m_sbpMode ?
"DG_INTEGRATE_GAUSS_LOBATTO" :
"DG_INTEGRATE_GAUS";
6433 dgIntegrationMethod =
6434 Context::getSolverProperty<MString>(
"dgIntegrationMethod", m_solverId, AT_, &dgIntegrationMethod);
6435 MString dgPolynomialType =
"DG_POLY_LEGENDRE";
6436 dgPolynomialType = Context::getSolverProperty<MString>(
"dgPolynomialType", m_solverId, AT_, &dgPolynomialType);
6439 parallelIo.setAttribute(dgIntegrationMethod,
"dgIntegrationMethod",
"polyDegs");
6440 parallelIo.setAttribute(dgPolynomialType,
"dgPolynomialType",
"polyDegs");
6443 for(
MInt i = 0; i < noVars; i++) {
6444 const MString name =
"variables" + to_string(i);
6445 parallelIo.defineArray(PIO_FLOAT, name, globalNoNodes);
6446 parallelIo.setAttribute(varNames[i],
"name", name);
6449 const MInt maxNoNodesXD =
ipow(m_maxNoNodes1D, nDim);
6450 const MInt elementDataSize = noVars * maxNoNodesXD;
6454 parallelIo.setOffset(grid().noInternalCells(), grid().domainOffset(domainId()));
6456 fill_n(&polyDegs[0], grid().noInternalCells(), m_initPolyDeg);
6457 for(
MInt elementId = 0; elementId < noElements; elementId++) {
6458 const MInt cellId = m_elements.cellId(elementId);
6459 polyDegs.
p[cellId] = m_elements.polyDeg(elementId);
6461 parallelIo.writeArray(polyDegs.
begin(),
"polyDegs");
6465 parallelIo.setOffset(grid().noInternalCells(), grid().domainOffset(domainId()));
6467 fill_n(&noNodes1D[0], grid().noInternalCells(), m_initNoNodes1D);
6468 for(
MInt elementId = 0; elementId < noElements; elementId++) {
6469 const MInt cellId = m_elements.cellId(elementId);
6470 noNodes1D.
p[cellId] = m_elements.noNodes1D(elementId);
6472 parallelIo.writeArray(noNodes1D.
begin(),
"noNodes1D");
6475 parallelIo.setOffset(localNoNodes, nodesOffset);
6476 for(
MInt i = 0; i < noVars; i++) {
6478 fill(buffer.
begin(), buffer.
end(), 0.0);
6479 for(
MInt e = 0; e < noElements; e++) {
6480 MFloat*
const b = &buffer[elementOffset[e]];
6481 const MInt dataOffset = e * elementDataSize + i;
6482 const MFloat*
const v = &data[dataOffset];
6483 for(
MInt n = 0; n < elementNodes[e]; n++) {
6484 b[n] = v[n * noVars];
6487 const MString name =
"variables" + to_string(i);
6488 parallelIo.writeArray(&buffer[0], name);
6500template <MInt nDim,
class SysEqn>
6505 m_log <<
"Saving file with node variable data for restarting... ";
6512 vector<MString> nodeVarNames(SysEqn::noNodeVars());
6513 for(
MInt nodeVar = 0; nodeVar < SysEqn::noNodeVars(); nodeVar++) {
6514 nodeVarNames[nodeVar] = m_sysEqn.nodeVarNames(nodeVar);
6518 saveNodalData(ss.str(), SysEqn::noNodeVars(), nodeVarNames, &m_elements.nodeVars(0));
6520 m_log <<
"done" << endl;
6531template <MInt nDim,
class SysEqn>
6536 m_log <<
"Loading nodevars file... ";
6540 ss << restartDir() <<
"nodevars" << getIdentifier(
g_multiSolverGrid,
"_b",
"") << ParallelIo::fileExt();
6543 using namespace parallel_io;
6544 ParallelIo parallelIo(ss.str(), PIO_READ, mpiComm());
6547 const MInt noElements = m_elements.size();
6549 for(
MInt cellId = 0, offset = 0, elementId = 0; elementId < noElements; cellId++) {
6552 if(cellId != m_elements.cellId(elementId)) {
6553 offset +=
ipow(m_initNoNodes1D, nDim);
6558 elementOffset[elementId] = offset;
6561 const MInt noNodesXD = m_elements.noNodesXD(elementId);
6562 offset += noNodesXD;
6572 const MInt localNoNodes = m_noTotalNodesXD + (grid().noInternalCells() - noElements) *
ipow(m_initNoNodes1D, nDim);
6575 ParallelIo::size_type nodesOffset, globalNoNodes;
6576 ParallelIo::calcOffset(localNoNodes, &nodesOffset, &globalNoNodes, mpiComm());
6579 parallelIo.setOffset(localNoNodes, nodesOffset);
6580 for(
MInt i = 0; i < m_sysEqn.noNodeVars(); i++) {
6582 const MString name =
"variables" + to_string(i);
6583 parallelIo.readArray(&buffer[0], name);
6584 for(
MInt e = 0; e < noElements; e++) {
6585 const MFloat*
const b = &buffer[elementOffset[e]];
6586 MFloat*
const v = &m_elements.nodeVars(e) + i;
6587 const MInt noNodesXD = m_elements.noNodesXD(e);
6588 for(
MInt n = 0; n < noNodesXD; n++) {
6589 v[n * m_sysEqn.noNodeVars()] =
b[n];
6594 m_log <<
"done" << endl;
6609template <MInt nDim,
class SysEqn>
6611 const MInt useNonSpecifiedRestartFile) {
6615 if(useNonSpecifiedRestartFile) {
6616 ss << restartDir() <<
"restart_" << getIdentifier(
g_multiSolverGrid,
"b",
"") << ParallelIo::fileExt();
6618 ss << restartDir() <<
"restart_" << getIdentifier(
g_multiSolverGrid,
"b") << setw(8) << setfill(
'0') << timeStep
6619 << ParallelIo::fileExt();
6621 const MString fileName = ss.str();
6627template <MInt nDim,
class SysEqn>
6631 BC bc = m_boundaryConditionFactory.create(bcId);
6646template <MInt nDim,
class SysEqn>
6649 RECORD_TIMER_START(m_timers[Timers::TimeDeriv]);
6659 prolongToSurfaces();
6662 applyForwardProjection();
6665 startMpiSurfaceExchange();
6669 calcVolumeIntegral();
6672 calcInnerSurfaceFlux();
6676 finishMpiSurfaceExchange();
6682 calcBoundarySurfaceFlux(tStage);
6685 calcMpiSurfaceFlux();
6688 calcSurfaceIntegral();
6694 calcSourceTerms(tStage);
6697 applyExternalSourceTerms(tStage);
6701 RECORD_TIMER_START(m_timers[Timers::Sponge]);
6702 sponge().calcSourceTerms();
6703 RECORD_TIMER_STOP(m_timers[Timers::Sponge]);
6706 RECORD_TIMER_STOP(m_timers[Timers::TimeDeriv]);
6713template <MInt nDim,
class SysEqn>
6716 RECORD_TIMER_START(m_timers[Timers::Prolong]);
6718 prolongToSurfaces(0, m_elements.size());
6720 RECORD_TIMER_STOP(m_timers[Timers::Prolong]);
6732template <MInt nDim,
class SysEqn>
6736 const MInt* polyDegs = &m_elements.polyDeg(0);
6737 const MInt* noNodes = &m_elements.noNodes1D(0);
6738 const MInt* surfaceIds = &m_elements.surfaceIds(0, 0);
6740 using namespace dg::interpolation;
6749#pragma omp parallel for
6752 for(
MInt elementId = begin; elementId < end; elementId++) {
6753 const MInt surfaceIdOffset = elementId * 2 * nDim;
6754 const MInt polyDeg = polyDegs[elementId];
6755 const MInt noNodes1D = noNodes[elementId];
6759 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
6760 const MInt srfcId = surfaceIds[surfaceIdOffset + dir];
6761 const MInt side = 1 - dir % 2;
6763 MFloat* src = &m_elements.variables(elementId);
6764 MFloat* dest = &m_surfaces.variables(srfcId, side);
6766 prolongToFaceGauss<nDim, SysEqn::noVars()>(src, dir, noNodes1D, &interp.
m_LFace[0][0], &interp.
m_LFace[1][0],
6772#pragma omp parallel for
6775 for(
MInt elementId = begin; elementId < end; elementId++) {
6776 const MInt surfaceIdOffset = elementId * 2 * nDim;
6777 const MInt noNodes1D = noNodes[elementId];
6780 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
6781 const MInt srfcId = surfaceIds[surfaceIdOffset + dir];
6782 const MInt side = 1 - dir % 2;
6784 MFloat* src = &m_elements.variables(elementId);
6785 MFloat* dest = &m_surfaces.variables(srfcId, side);
6787 prolongToFaceGaussLobatto<nDim, SysEqn::noVars()>(src, dir, noNodes1D, dest);
6804template <MInt nDim,
class SysEqn>
6807 RECORD_TIMER_START(m_timers[Timers::ForwardProjection]);
6814 const MInt* surfaceIds = &m_elements.surfaceIds(0, 0);
6815 const MInt noVars = SysEqn::noVars();
6816 const MInt maxNoNodes1D = m_maxNoNodes1D;
6817 const MInt maxNoNodes1D3 = (nDim == 3) ? maxNoNodes1D : 1;
6818 const MInt noHElements = m_helements.size();
6819 const MInt noDirs = 2 * nDim;
6820 const MInt noSurfs = 2 * (nDim - 1);
6826 if(noHElements > 0) {
6828#pragma omp parallel for
6830 for(
MInt hElementId = 0; hElementId < noHElements; hElementId++) {
6831 const MInt coarseElementId = m_helements.elementId(hElementId);
6832 for(
MInt dir = 0; dir < noDirs; dir++) {
6833 const MInt coarseSrfcId = m_elements.surfaceIds(coarseElementId, dir);
6834 for(
MInt pos = 0; pos < noSurfs; pos++) {
6835 const MInt hSrfcId = m_helements.hrefSurfaceIds(hElementId, dir, pos);
6836 const MInt side = 1 - dir % 2;
6844 if(coarseSrfcId != hSrfcId) {
6845 const MInt noNodes1D = m_elements.noNodes1D(coarseElementId);
6846 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
6847 const MInt size = noNodes1D * noNodes1D3 * SysEqn::noVars();
6850 copy_n(&m_surfaces.variables(coarseSrfcId, side), size, &m_surfaces.variables(hSrfcId, side));
6860 const MInt noElements = m_elements.size();
6861 MFloatTensor projected(maxNoNodes1D, maxNoNodes1D3, noVars);
6864 if(noHElements > 0) {
6866#pragma omp parallel for firstprivate(projected)
6868 for(
MInt hElementId = 0; hElementId < noHElements; hElementId++) {
6869 for(
MInt dir = 0; dir < noDirs; dir++) {
6870 for(
MInt pos = 0; pos < noSurfs; pos++) {
6871 const MInt hSrfcId = m_helements.hrefSurfaceIds(hElementId, dir, pos);
6872 const MInt side = 1 - dir % 2;
6879 const MInt surfaceNoNodes1D = m_surfaces.noNodes1D(hSrfcId);
6880 const MInt surfaceNoNodes1D3 = (nDim == 3) ? surfaceNoNodes1D : 1;
6881 const MInt size = surfaceNoNodes1D * surfaceNoNodes1D3 * noVars;
6884 calcMortarProjection<dg::mortar::forward, SysEqn::noVars()>(
6885 hSrfcId, dir, &m_surfaces.variables(hSrfcId, side), &projected[0], m_elements, m_surfaces);
6888 copy_n(&projected[0], size, &m_surfaces.variables(hSrfcId, side));
6896#pragma omp parallel for firstprivate(projected)
6898 for(
MInt elementId = 0; elementId < noElements; elementId++) {
6899 const MInt surfaceIdOffset = elementId * 2 * nDim;
6901 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
6903 const MInt srfcId = surfaceIds[surfaceIdOffset + dir];
6906 if(srfcId < m_innerSurfacesOffset) {
6910 const MInt surfacePolyDeg = m_surfaces.polyDeg(srfcId);
6911 const MInt elementPolyDeg = m_elements.polyDeg(elementId);
6912 const MInt side = 1 - dir % 2;
6913 const MInt surfaceNoNodes1D = m_surfaces.noNodes1D(srfcId);
6914 const MInt surfaceNoNodes1D3 = (nDim == 3) ? surfaceNoNodes1D : 1;
6915 const MInt size = surfaceNoNodes1D * surfaceNoNodes1D3 * noVars;
6918 if(m_surfaces.fineCellId(srfcId) != -1 && m_surfaces.fineCellId(srfcId) != m_elements.cellId(elementId)) {
6923 if(surfacePolyDeg > elementPolyDeg) {
6925 calcMortarProjection<dg::mortar::forward, SysEqn::noVars()>(srfcId, dir, &m_surfaces.variables(srfcId, side),
6926 &projected[0], m_elements, m_surfaces);
6929 copy_n(&projected[0], size, &m_surfaces.variables(srfcId, side));
6934 RECORD_TIMER_STOP(m_timers[Timers::ForwardProjection]);
6943template <MInt nDim,
class SysEqn>
6946 RECORD_TIMER_START(m_timers[Timers::VolInt]);
6948 calcVolumeIntegral(m_elements.size(), m_elements, m_sysEqn);
6950 RECORD_TIMER_STOP(m_timers[Timers::VolInt]);
6964template <MInt nDim,
class SysEqn>
6967 RECORD_TIMER_START(m_timers[Timers::Flux]);
6968 RECORD_TIMER_START(m_timers[Timers::FluxBndry]);
6970 for(
auto&& bc : m_boundaryConditions) {
6974 RECORD_TIMER_STOP(m_timers[Timers::FluxBndry]);
6975 RECORD_TIMER_STOP(m_timers[Timers::Flux]);
6986template <MInt nDim,
class SysEqn>
6989 RECORD_TIMER_START(m_timers[Timers::Flux]);
6990 RECORD_TIMER_START(m_timers[Timers::FluxInner]);
6992 const MInt begin = m_innerSurfacesOffset;
6993 const MInt end = m_innerSurfacesOffset + m_noInnerSurfaces;
6994 calcRegularSurfaceFlux(begin, end, m_surfaces, m_sysEqn);
6996 RECORD_TIMER_STOP(m_timers[Timers::FluxInner]);
6997 RECORD_TIMER_STOP(m_timers[Timers::Flux]);
7008template <MInt nDim,
class SysEqn>
7011 RECORD_TIMER_START(m_timers[Timers::Flux]);
7012 RECORD_TIMER_START(m_timers[Timers::FluxMPI]);
7014 const MInt begin = m_mpiSurfacesOffset;
7015 const MInt end = m_mpiSurfacesOffset + m_noMpiSurfaces;
7016 calcRegularSurfaceFlux(begin, end, m_surfaces, m_sysEqn);
7018 RECORD_TIMER_STOP(m_timers[Timers::FluxMPI]);
7019 RECORD_TIMER_STOP(m_timers[Timers::Flux]);
7028template <MInt nDim,
class SysEqn>
7031 RECORD_TIMER_START(m_timers[Timers::SurfInt]);
7033 calcSurfaceIntegral(0, m_elements.size(), m_elements, m_surfaces, m_helements, m_helements.size());
7035 RECORD_TIMER_STOP(m_timers[Timers::SurfInt]);
7050template <MInt nDim,
class SysEqn>
7056 const MInt noHElements) {
7059 const MInt noVars = SysEqn::noVars();
7068#pragma omp parallel for
7071 for(
MInt elementId = begin; elementId < end; elementId++) {
7072 const MInt surfaceIdOffset = elementId * 2 * nDim;
7074 for(
MInt dir = 0; dir < 2 * nDim; dir++) {
7076 const MInt srfcId = surfaceIds[surfaceIdOffset + dir];
7079 const MInt srfcNodes1D3 = (nDim == 3) ? srfcNodes1D : 1;
7084 const MInt side = 1 - dir % 2;
7092 if(srfcPolyDeg > polyDeg) {
7096 calcMortarProjection<dg::mortar::reverse, SysEqn::noVars()>(srfcId, dir, &f[0], &fp[0], elem, surf);
7099 applySurfaceIntegral(&elem.
rightHandSide(elementId), polyDeg, noNodes1D, srfcId, side, &fp[0], surf);
7103 applySurfaceIntegral(&elem.
rightHandSide(elementId), polyDeg, noNodes1D, srfcId, side, &f[0], surf);
7111 if(noHElements > 0) {
7112 const MInt noDirs = 2 * nDim;
7113 const MInt noSurfs = 2 * (nDim - 1);
7116#pragma omp parallel for
7118 for(
MInt hElementId = 0; hElementId < noHElements; hElementId++) {
7120 const MInt polyDeg = elem.
polyDeg(coarseElementId);
7123 for(
MInt dir = 0; dir < noDirs; dir++) {
7126 for(
MInt pos = 0; pos < noSurfs; pos++) {
7128 const MInt side = 1 - dir % 2;
7136 const MInt srfcNodes1D3 = (nDim == 3) ? srfcNodes1D : 1;
7141 calcMortarProjection<dg::mortar::reverse, SysEqn::noVars()>(hSrfcId, dir, &f[0], &fp[0], elem, surf);
7144 applySurfaceIntegral(&elem.
rightHandSide(coarseElementId), polyDeg, noNodes1D, coarseSrfcId, side, &fp[0],
7168template <MInt nDim,
class SysEqn>
7173 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
7174 const MInt noVars = SysEqn::noVars();
7176 const MInt srfcNodes1D3 = (nDim == 3) ? srfcNodes1D : 1;
7181 MFloatTensor ut(rhs, noNodes1D, noNodes1D, noNodes1D3, noVars);
7183 const MInt index = (side == 0) ? (noNodes1D - 1) : 0;
7184 const MInt sign = (side == 0) ? 1 : -1;
7192 for(
MInt i = 0; i < noNodes1D; i++) {
7193 for(
MInt j = 0; j < noNodes1D; j++) {
7194 for(
MInt k = 0; k < noNodes1D3; k++) {
7195 for(
MInt n = 0; n < noVars; n++) {
7196 ut(i, j, k, n) += sign * f(j, k, n) * interp.
m_LhatFace[1 - side][i];
7204 for(
MInt i = 0; i < noNodes1D; i++) {
7205 for(
MInt j = 0; j < noNodes1D; j++) {
7206 for(
MInt k = 0; k < noNodes1D3; k++) {
7207 for(
MInt n = 0; n < noVars; n++) {
7208 ut(i, j, k, n) += sign * f(i, k, n) * interp.
m_LhatFace[1 - side][j];
7216 for(
MInt i = 0; i < noNodes1D; i++) {
7217 for(
MInt j = 0; j < noNodes1D; j++) {
7218 for(
MInt k = 0; k < noNodes1D3; k++) {
7219 for(
MInt n = 0; n < noVars; n++) {
7220 ut(i, j, k, n) += sign * f(i, j, n) * interp.
m_LhatFace[1 - side][k];
7227 mTerm(1, AT_,
"Bad direction id");
7233 for(
MInt j = 0; j < noNodes1D; j++) {
7234 for(
MInt k = 0; k < noNodes1D3; k++) {
7235 for(
MInt n = 0; n < noVars; n++) {
7236 ut(index, j, k, n) += sign * f(j, k, n) * interp.
m_LhatFace[1 - side][index];
7243 for(
MInt i = 0; i < noNodes1D; i++) {
7244 for(
MInt k = 0; k < noNodes1D3; k++) {
7245 for(
MInt n = 0; n < noVars; n++) {
7246 ut(i, index, k, n) += sign * f(i, k, n) * interp.
m_LhatFace[1 - side][index];
7253 for(
MInt i = 0; i < noNodes1D; i++) {
7254 for(
MInt j = 0; j < noNodes1D3; j++) {
7255 for(
MInt n = 0; n < noVars; n++) {
7256 ut(i, j, index, n) += sign * f(i, j, n) * interp.
m_LhatFace[1 - side][index];
7263 mTerm(1, AT_,
"Bad direction id");
7273template <MInt nDim,
class SysEqn>
7276 RECORD_TIMER_START(m_timers[Timers::Jacobian]);
7278 applyJacobian(m_elements.size(), m_elements);
7280 RECORD_TIMER_STOP(m_timers[Timers::Jacobian]);
7293template <MInt nDim,
class SysEqn>
7300#pragma omp parallel for
7302 for(
MInt elementId = 0; elementId < noElements; elementId++) {
7303 const MInt dataBlockSize = elem.
noNodesXD(elementId) * SysEqn::noVars();
7304 const MFloat invJacobian = invJacobians[elementId];
7307 for(
MInt dataId = 0; dataId < dataBlockSize; dataId++) {
7308 rhs[dataId] *= -invJacobian;
7319template <MInt nDim,
class SysEqn>
7322 RECORD_TIMER_START(m_timers[Timers::Sources]);
7324 calcSourceTerms(t, m_elements.size(), m_elements, m_sysEqn);
7326 RECORD_TIMER_STOP(m_timers[Timers::Sources]);
7331template <MInt nDim,
class SysEqn>
7334 RECORD_TIMER_START(m_timers[Timers::ExternalSources]);
7336 const MInt* noNodes1D = &m_elements.noNodes1D(0);
7339 const MFloat rampUpFactor =
7343#pragma omp parallel for
7345 for(
MInt elementId = 0; elementId < m_elements.size(); elementId++) {
7346 const MInt dataBlockSize =
ipow(noNodes1D[elementId], nDim) * SysEqn::noVars();
7347 MFloat*
const rhs = &m_elements.rightHandSide(elementId);
7348 const MFloat*
const sources = &m_elements.externalSource(elementId);
7350 for(
MInt dataId = 0; dataId < dataBlockSize; dataId++) {
7351 rhs[dataId] += rampUpFactor * sources[dataId];
7355 RECORD_TIMER_STOP(m_timers[Timers::ExternalSources]);
7365template <MInt nDim,
class SysEqn>
7367 return (domainId() == 0);
7378template <MInt nDim,
class SysEqn>
7387 RECORD_TIMER_START(m_timers[Timers::SurfExchange]);
7388 RECORD_TIMER_START(m_timers[Timers::Accumulated]);
7389 RECORD_TIMER_START(m_timers[Timers::MPI]);
7392 RECORD_TIMER_START(m_timers[Timers::MPIComm]);
7393 RECORD_TIMER_START(m_timers[Timers::SurfExchangeComm]);
7394 RECORD_TIMER_START(m_timers[Timers::SECommRecv]);
7395 if(!m_mpiRecvRequestsOpen && m_noExchangeNghbrDomains > 0) {
7396 MPI_Startall(m_noExchangeNghbrDomains, &m_recvRequests[0], AT_);
7397 m_mpiRecvRequestsOpen =
true;
7399 RECORD_TIMER_STOP(m_timers[Timers::SECommRecv]);
7400 RECORD_TIMER_STOP(m_timers[Timers::SurfExchangeComm]);
7401 RECORD_TIMER_STOP(m_timers[Timers::MPIComm]);
7404 RECORD_TIMER_START(m_timers[Timers::MPIWait]);
7405 RECORD_TIMER_START(m_timers[Timers::SurfExchangeWait]);
7406 RECORD_TIMER_START(m_timers[Timers::SEWaitSend]);
7407 if(m_mpiSendRequestsOpen && m_noExchangeNghbrDomains > 0) {
7408 MPI_Waitall(m_noExchangeNghbrDomains, &m_sendRequests[0], MPI_STATUSES_IGNORE, AT_);
7409 m_mpiSendRequestsOpen =
false;
7411 RECORD_TIMER_STOP(m_timers[Timers::SEWaitSend]);
7412 RECORD_TIMER_STOP(m_timers[Timers::SurfExchangeWait]);
7413 RECORD_TIMER_STOP(m_timers[Timers::MPIWait]);
7416#ifdef DG_USE_MPI_BUFFERS
7417 RECORD_TIMER_START(m_timers[Timers::MPICopy]);
7418 RECORD_TIMER_START(m_timers[Timers::SurfExchangeCopy]);
7419 RECORD_TIMER_START(m_timers[Timers::SECopySend]);
7422 const MInt dataBlockSize =
ipow(m_maxNoNodes1D, nDim - 1) * SysEqn::noVars();
7424 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
7426 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
7427 const MInt srfcId = m_mpiSurfaces[i][j];
7428 const MInt sideId = m_surfaces.internalSideId(srfcId);
7431 const MFloat* data = &m_surfaces.variables(srfcId, sideId);
7432 copy(data, data + dataBlockSize, &m_sendBuffers[i][size]);
7433 size += dataBlockSize;
7435 ASSERT(size ==
static_cast<MInt>(m_sendBuffers[i].size()),
"Data size does not match buffer size.");
7437 RECORD_TIMER_STOP(m_timers[Timers::SECopySend]);
7438 RECORD_TIMER_STOP(m_timers[Timers::SurfExchangeCopy]);
7439 RECORD_TIMER_STOP(m_timers[Timers::MPICopy]);
7443 RECORD_TIMER_START(m_timers[Timers::MPIComm]);
7444 RECORD_TIMER_START(m_timers[Timers::SurfExchangeComm]);
7445 RECORD_TIMER_START(m_timers[Timers::SECommSend]);
7446 if(m_noExchangeNghbrDomains > 0) {
7447 MPI_Startall(m_noExchangeNghbrDomains, &m_sendRequests[0], AT_);
7449 RECORD_TIMER_STOP(m_timers[Timers::SECommSend]);
7450 RECORD_TIMER_STOP(m_timers[Timers::SurfExchangeComm]);
7451 RECORD_TIMER_STOP(m_timers[Timers::MPIComm]);
7453 m_mpiSendRequestsOpen =
true;
7455 RECORD_TIMER_STOP(m_timers[Timers::MPI]);
7456 RECORD_TIMER_STOP(m_timers[Timers::Accumulated]);
7457 RECORD_TIMER_STOP(m_timers[Timers::SurfExchange]);
7468template <MInt nDim,
class SysEqn>
7471 RECORD_TIMER_START(m_timers[Timers::SurfExchange]);
7472 RECORD_TIMER_START(m_timers[Timers::Accumulated]);
7473 RECORD_TIMER_START(m_timers[Timers::MPI]);
7480 if(noDomains() > 1 && (!m_mpiRecvRequestsOpen || !m_mpiSendRequestsOpen)) {
7481 TERMM(1,
"MPI requests are not open: receive=" + std::to_string(m_mpiRecvRequestsOpen)
7482 +
", send=" + std::to_string(m_mpiSendRequestsOpen));
7485#ifdef DG_USE_MPI_WAITSOME
7486 TERMM(1,
"This code is untested and will not work with p-refinement.");
7489 vector<MInt> indices(m_noExchangeNghbrDomains);
7494 RECORD_TIMER_START(m_timers[Timers::MPIWait]);
7495 RECORD_TIMER_START(m_timers[Timers::SurfExchangeWait]);
7496 RECORD_TIMER_START(m_timers[Timers::SEWaitRecv]);
7497 MPI_Waitsome(m_noExchangeNghbrDomains, &m_recvRequests[0], &noFinished, &indices[0], MPI_STATUSES_IGNORE, AT_);
7498 while(noFinished != MPI_UNDEFINED) {
7499 for(
MInt i = 0; i < noFinished; i++) {
7500 const MInt index = indices[i];
7502 for(
size_t j = 0; j < m_mpiSurfaces[index].size(); j++) {
7503 const MInt srfcId = m_mpiSurfaces[index][j];
7504 const MInt noNodes1D = m_surfaces.noNodes1D(srfcId);
7505 const MInt noNodesXD =
ipow(noNodes1D, nDim - 1);
7506 const MInt sideId = 1 - m_surfaces.internalSideId(srfcId);
7508 MFloat* data = m_surfaces.variables(srfcId, sideId);
7509 const MInt dataBlockSize = noNodesXD * SysEqn::noVars();
7510 copy_n(&m_recvBuffers[index][size], dataBlockSize, data);
7511 size += dataBlockSize;
7513 ASSERT(size ==
static_cast<MInt>(m_recvBuffers[index].size()),
"Data size does not match buffer size.");
7517 MPI_Waitsome(m_noExchangeNghbrDomains, &m_recvRequests[0], &noFinished, &indices[0], MPI_STATUSES_IGNORE, AT_);
7519 RECORD_TIMER_STOP(m_timers[Timers::SEWaitRecv]);
7520 RECORD_TIMER_STOP(m_timers[Timers::SurfExchangeWait]);
7521 RECORD_TIMER_STOP(m_timers[Timers::MPIWait]);
7525 RECORD_TIMER_START(m_timers[Timers::MPIWait]);
7526 RECORD_TIMER_START(m_timers[Timers::SurfExchangeWait]);
7527 RECORD_TIMER_START(m_timers[Timers::SEWaitRecv]);
7528 if(m_noExchangeNghbrDomains > 0) {
7529 MPI_Waitall(m_noExchangeNghbrDomains, &m_recvRequests[0], MPI_STATUSES_IGNORE, AT_);
7531 RECORD_TIMER_STOP(m_timers[Timers::SEWaitRecv]);
7532 RECORD_TIMER_STOP(m_timers[Timers::SurfExchangeWait]);
7533 RECORD_TIMER_STOP(m_timers[Timers::MPIWait]);
7536#ifdef DG_USE_MPI_BUFFERS
7537 RECORD_TIMER_START(m_timers[Timers::MPICopy]);
7538 RECORD_TIMER_START(m_timers[Timers::SurfExchangeCopy]);
7539 RECORD_TIMER_START(m_timers[Timers::SECopyRecv]);
7542 const MInt dataBlockSize =
ipow(m_maxNoNodes1D, nDim - 1) * SysEqn::noVars();
7544 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
7546 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
7549 const MInt srfcId = m_mpiRecvSurfaces[i][j];
7553 const MInt sideId = 1 - m_surfaces.internalSideId(srfcId);
7554 MFloat* data = &m_surfaces.variables(srfcId, sideId);
7555 copy_n(&m_recvBuffers[i][size], dataBlockSize, data);
7556 size += dataBlockSize;
7558 ASSERT(size ==
static_cast<MInt>(m_recvBuffers[i].size()),
"Data size does not match buffer size.");
7560 RECORD_TIMER_STOP(m_timers[Timers::SECopyRecv]);
7561 RECORD_TIMER_STOP(m_timers[Timers::SurfExchangeCopy]);
7562 RECORD_TIMER_STOP(m_timers[Timers::MPICopy]);
7566 m_mpiRecvRequestsOpen =
false;
7569 if(m_noExchangeNghbrDomains > 0) {
7570 MPI_Startall(m_noExchangeNghbrDomains, &m_recvRequests[0], AT_);
7571 m_mpiRecvRequestsOpen =
true;
7574 RECORD_TIMER_STOP(m_timers[Timers::MPI]);
7575 RECORD_TIMER_STOP(m_timers[Timers::Accumulated]);
7576 RECORD_TIMER_STOP(m_timers[Timers::SurfExchange]);
7586template <MInt nDim,
class SysEqn>
7589 RECORD_TIMER_START(m_timers[Timers::CalcTimeStep]);
7591 MFloat minDt = numeric_limits<MFloat>::infinity();
7594 const MBool forcedTimeStep = (m_externalDt > 0.0);
7599 const MInt noElements = m_elements.size();
7600 const MInt* noNodes1D = &m_elements.noNodes1D(0);
7601 const MFloat* invJacobians = &m_elements.invJacobian(0);
7604#pragma omp parallel for reduction(min : minDt)
7606 for(
MInt elementId = 0; elementId < noElements; elementId++) {
7607 MFloat ts = m_sysEqn.getTimeStep(&m_elements.nodeVars(elementId), &m_elements.variables(elementId),
7608 noNodes1D[elementId], invJacobians[elementId], m_sbpMode);
7609 minDt = min(minDt, ts);
7613 if(forcedTimeStep) {
7614 if(m_externalDt > minDt) {
7615 cerr <<
" WARNING: external dt larger than computed time step size on global domain " <<
globalDomainId() <<
"("
7616 << m_externalDt <<
" > " << minDt <<
")" << std::endl;
7621 m_log <<
"DG-Solver #" << solverId() <<
" using external time step size: " << m_externalDt
7622 <<
" (computed global minimum: " << minDt <<
")" << std::endl;
7623 minDt = m_externalDt;
7627 if(std::isnan(minDt)) {
7629 "Time step is NaN at time step " + to_string(m_timeStep) +
" on global domain "
7636 "Time step is less than zero at time step " + to_string(m_timeStep) +
" on global domain "
7642 if(minDt < 1.0e-100) {
7644 "Time step is less than 1.0e-100 at time step " + to_string(m_timeStep) +
" on global domain "
7652 RECORD_TIMER_STOP(m_timers[Timers::CalcTimeStep]);
7674template <MInt nDim,
class SysEqn>
7678 RECORD_TIMER_START(m_timers[Timers::RungeKuttaStep]);
7681 const MInt totalSize = m_internalDataSize;
7684 for(
MInt s = 0; s < m_noRkStages; s++) {
7686 if(substep != -1 && s != substep) {
7690 calcDgTimeDerivative(t, t + m_timeIntegrationCoefficientsC[s] * dt);
7692 RECORD_TIMER_START(m_timers[Timers::TimeInt]);
7693 subTimeStepRk(dt, s, totalSize, &m_elements.rightHandSide(0), &m_elements.variables(0),
7694 &m_elements.timeIntStorage(0));
7695 RECORD_TIMER_STOP(m_timers[Timers::TimeInt]);
7702 const MBool lastRkStage = (s == (m_noRkStages - 1));
7703 const MBool adaptationTimeStep = (lastRkStage && isAdaptationTimeStep(m_timeStep + 1));
7708 RECORD_TIMER_START(m_timers[Timers::TimeDeriv]);
7709 prolongToSurfaces();
7710 applyForwardProjection();
7711 startMpiSurfaceExchange();
7712 RECORD_TIMER_STOP(m_timers[Timers::TimeDeriv]);
7717 RECORD_TIMER_STOP(m_timers[Timers::RungeKuttaStep]);
7732template <MInt nDim,
class SysEqn>
7735 const MInt totalSize,
7738 MFloat*
const timeIntStorage) {
7742 MFloat*
const p = variables;
7743 MFloat*
const k = timeIntStorage;
7747 bDt = m_timeIntegrationCoefficientsB[stage] * dt;
7752#pragma omp parallel for
7754 for(
MInt i = 0; i < totalSize; i++) {
7761#pragma omp parallel for
7763 for(
MInt i = 0; i < totalSize; i++) {
7764 k[i] = rhs[i] - k[i] * m_timeIntegrationCoefficientsA[stage];
7780template <MInt nDim,
class SysEqn>
7784 RECORD_TIMER_START(m_timers[Timers::Accumulated]);
7785 RECORD_TIMER_START(m_timers[Timers::AnalyzeTimeStep]);
7787 vector<MFloat> L2Error(m_sysEqn.noVars());
7788 vector<MFloat> LInfError(m_sysEqn.noVars());
7789 vector<MFloat> L2ErrLocal(m_sysEqn.noVars());
7790 vector<MFloat> LInfErrLocal(m_sysEqn.noVars());
7791 if(m_calcErrorNorms) {
7792 calcErrorNorms(t, L2Error, LInfError, L2ErrLocal, LInfErrLocal);
7795 if(any_of(L2Error.begin(), L2Error.end(), [](
const MFloat e) { return std::isnan(e); })
7796 || any_of(LInfError.begin(), LInfError.end(), [](
const MFloat e) { return std::isnan(e); })) {
7797 saveSolutionFile(
"nan_" + to_string(timeStep));
7801 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
7802 if(std::isnan(L2ErrLocal[i])) {
7803 TERMM(1,
"L^2 error for variable '" + m_sysEqn.consVarNames(i) +
"' is NaN at time step " + to_string(timeStep)
7806 if(std::isnan(LInfErrLocal[i])) {
7807 TERMM(1,
"L^inf error for variable '" + m_sysEqn.consVarNames(i) +
"' is NaN at time step "
7808 + to_string(timeStep) +
" on global domain " + to_string(
globalDomainId()));
7813 const MInt precision = m_noErrorDigits;
7814 const MInt fieldwith = precision + 6;
7820 log <<
"----------------------------------------"
7821 <<
"----------------------------------------" << endl;
7822 log <<
" Solver " << solverId() <<
" running '" << m_sysEqn.sysEqnName() <<
"' with N = " << m_initPolyDeg
7823 <<
" and maxLevel = " << grid().maxLevel() << endl;
7824 log <<
"----------------------------------------"
7825 <<
"----------------------------------------" << endl;
7826 log <<
" No. timesteps: " << timeStep << endl;
7827 log <<
" dt: " << scientific << dt << endl;
7828 log <<
" Run time: " << runTimeTotal <<
" s" << endl;
7829 log <<
" Time/DOF/step: " << runTimeRelative <<
" s" << endl;
7831 if(m_calcErrorNorms) {
7832 streamsize ss = log.precision();
7833 log << setprecision(precision);
7834 log <<
" Variable: ";
7835 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
7836 log <<
" " << left << setw(fieldwith) << setfill(
' ') << m_sysEqn.consVarNames(i);
7838 log << right << endl;
7839 log <<
" L^2 error: ";
7840 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
7841 log <<
" " << L2Error[i];
7844 log <<
" L^inf error: ";
7845 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
7846 log <<
" " << LInfError[i];
7849 log << setprecision(ss);
7852 log <<
"----------------------------------------"
7853 <<
"----------------------------------------" << endl;
7854 log <<
" Simulation time: " << t << endl;
7855 log <<
"----------------------------------------"
7856 <<
"----------------------------------------" << endl;
7858 m_log << log.str() << endl;
7859 cout << log.str() << endl;
7862 MBool finalTimeStep =
false;
7863 if(m_finalTime - m_time - m_dt < 1.0E-10) {
7864 finalTimeStep =
true;
7865 }
else if(m_timeStep == m_timeSteps) {
7866 finalTimeStep =
true;
7869 if(m_calcErrorNorms && isMpiRoot() && finalTimeStep) {
7870 std::ofstream eocStream;
7871 eocStream.
open(
"EOCout.csv");
7872 eocStream << m_statGlobalNoActiveDOFs <<
"," << L2Error[0] <<
"," << LInfError[0];
7875 std::ofstream perfStream;
7876 perfStream.open(
"TimerOut.csv");
7877 perfStream << m_statGlobalNoActiveCells <<
"," << m_timeStep <<
"," << runTimeTotal <<
"," << runTimeRelative;
7883 if(m_calcErrorNorms) {
7884 stringstream logLocal;
7885 logLocal << setprecision(precision) << scientific;
7886 logLocal << endl <<
" Variable: ";
7887 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
7888 logLocal <<
" " << left << setw(fieldwith) << setfill(
' ') << m_sysEqn.consVarNames(i);
7890 logLocal << right << endl;
7891 logLocal <<
" Local L^2 error: ";
7892 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
7893 logLocal <<
" " << L2ErrLocal[i];
7896 logLocal <<
" Local L^inf error:";
7897 for(
MInt i = 0; i < m_sysEqn.noVars(); i++) {
7898 logLocal <<
" " << LInfErrLocal[i];
7902 m_log << logLocal.str() << endl;
7905 RECORD_TIMER_STOP(m_timers[Timers::AnalyzeTimeStep]);
7906 RECORD_TIMER_STOP(m_timers[Timers::Accumulated]);
7924template <MInt nDim,
class SysEqn>
7926 vector<MFloat>& L2Error,
7927 vector<MFloat>& LInfError,
7928 vector<MFloat>& L2ErrLocal,
7929 vector<MFloat>& LInfErrLocal) {
7931 const MInt noVars = m_sysEqn.noVars();
7932 L2Error.assign(noVars, F0);
7933 LInfError.assign(noVars, F0);
7934 L2ErrLocal.assign(noVars, F0);
7935 LInfErrLocal.assign(noVars, F0);
7939 const MInt noNodesAnalysis1D = m_noAnalysisNodes;
7940 const MInt noNodesAnalysis1D3 = (nDim == 3) ? noNodesAnalysis1D : 1;
7941 MFloatTensor u(noNodesAnalysis1D, noNodesAnalysis1D, noNodesAnalysis1D3, noVars);
7942 MFloatTensor x(noNodesAnalysis1D, noNodesAnalysis1D, noNodesAnalysis1D3, nDim);
7945 MFloatTensor nodeVars(noNodesAnalysis1D, noNodesAnalysis1D, noNodesAnalysis1D3, max(SysEqn::noNodeVars(), 1));
7948 const MInt noElements = m_elements.size();
7949 for(
MInt elementId = 0; elementId < noElements; elementId++) {
7951 const MInt polyDegElement = m_elements.polyDeg(elementId);
7952 const MInt noNodesElement = m_elements.noNodes1D(elementId);
7953 auto vdm = m_vdmAnalysis[polyDegElement][noNodesElement];
7954 dg::interpolation::interpolateNodes<nDim>(&m_elements.variables(elementId), &vdm[0], noNodesElement,
7955 m_noAnalysisNodes, noVars, &u[0]);
7958 dg::interpolation::interpolateNodes<nDim>(&m_elements.nodeCoordinates(elementId), &vdm[0], noNodesElement,
7959 m_noAnalysisNodes, nDim, &x[0]);
7962 const MFloat jacobian = pow(F1 / m_elements.invJacobian(elementId), nDim);
7963 for(
MInt i = 0; i < noNodesAnalysis1D; i++) {
7964 for(
MInt j = 0; j < noNodesAnalysis1D; j++) {
7965 for(
MInt k = 0; k < noNodesAnalysis1D3; k++) {
7968 m_sysEqn.calcInitialCondition(t, &x(i, j, k, 0), &nodeVars(i, j, k, 0), &uExact[0]);
7969 for(
MInt v = 0; v < noVars; v++) {
7970 const MFloat diff = uExact[v] - u(i, j, k, v);
7971 L2ErrLocal[v] += pow(diff, F2) * m_wVolumeAnalysis(i, j, k) * jacobian;
7972 LInfErrLocal[v] = max(LInfErrLocal[v], fabs(diff));
7980 RECORD_TIMER_START(m_timers[Timers::MPI]);
7981 RECORD_TIMER_START(m_timers[Timers::MPIComm]);
7982 MPI_Allreduce(&L2ErrLocal[0], &L2Error[0], noVars, MPI_DOUBLE, MPI_SUM, mpiComm(), AT_,
"L2ErrLocal[0]",
7984 MPI_Allreduce(&LInfErrLocal[0], &LInfError[0], noVars, MPI_DOUBLE, MPI_MAX, mpiComm(), AT_,
"LInfErrLocal[0]",
7986 RECORD_TIMER_STOP(m_timers[Timers::MPIComm]);
7987 RECORD_TIMER_STOP(m_timers[Timers::MPI]);
7990 for(
MInt v = 0; v < noVars; v++) {
7991 L2ErrLocal[v] = sqrt(L2ErrLocal[v] / m_localVolume);
7992 L2Error[v] = sqrt(L2Error[v] / m_globalVolume);
8001template <MInt nDim,
class SysEqn>
8004 RECORD_TIMER_START(m_timers[Timers::ResetRHS]);
8006 resetBuffer(m_internalDataSize, &m_elements.rightHandSide(0));
8008 RECORD_TIMER_STOP(m_timers[Timers::ResetRHS]);
8012template <MInt nDim,
class SysEqn>
8015 RECORD_TIMER_START(m_timers[Timers::ResetExternalSources]);
8017 resetBuffer(m_internalDataSize, &m_elements.externalSource(0));
8019 RECORD_TIMER_STOP(m_timers[Timers::ResetExternalSources]);
8032template <MInt nDim,
class SysEqn>
8037#pragma omp parallel for
8039 for(
MInt i = 0; i < totalSize; i++) {
8060template <MInt nDim,
class SysEqn>
8064 MInt foundElementId = -1;
8065 const MInt noElements = m_elements.size();
8067 for(
MInt elementId = 0; elementId < noElements; elementId++) {
8068 const MInt cellId = m_elements.cellId(elementId);
8069 const MFloat cellLength = grid().cellLengthAtCell(cellId);
8070 MBool pointInCell =
true;
8071 MBool takeCell =
true;
8073 for(
MInt i = 0; i < nDim; i++) {
8074 if(!(fabs(grid().tree().coordinate(cellId, i) - point[i]) <= cellLength * 0.5)) {
8075 pointInCell =
false;
8082 if(pointInCell && globalUnique) {
8083 for(
MInt dimId = 0; dimId < nDim; dimId++) {
8084 const MFloat distance = grid().tree().coordinate(cellId, dimId) - point[dimId];
8086 if(
approx(fabs(distance), cellLength * 0.5, MFloatEps)) {
8088 if(distance > 0.0) {
8091 const MInt neighborDir = 2 * dimId;
8092 const MInt surfaceId = m_elements.surfaceIds(elementId, neighborDir);
8095 if(isMpiSurface(surfaceId)) {
8103 if(pointInCell && takeCell) {
8104 foundElementId = elementId;
8109 return foundElementId;
8142template <MInt nDim,
class SysEqn>
8145 const MInt elementId = getElementIdAtPoint(point);
8147 if(elementId == -1) {
8151 calcStateAtPoint(point, elementId, state);
8170template <MInt nDim,
class SysEqn>
8175 const MInt cellId = m_elements.cellId(elementId);
8176 const MInt polyDeg = m_elements.polyDeg(elementId);
8177 const MInt noNodes1D = m_elements.noNodes1D(elementId);
8178 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
8181 const MFloat cellLength = grid().cellLengthAtCell(cellId);
8184 MFloat pointOnUnitInt[nDim];
8186 const MFloatTensor U(
const_cast<MFloat*
>(u), noNodes1D, noNodes1D, noNodes1D3, noVars);
8190 for(
MInt n = 0; n < nDim; n++) {
8191 pointOnUnitInt[n] = F2 * (point[n] - grid().tree().coordinate(cellId, n)) / cellLength;
8195 IF_CONSTEXPR(nDim == 2) {
8212 for(
MInt n = 0; n < nDim; n++) {
8213 lagrangePoly[n].
resize(noNodes1D);
8214 fill_n(&lagrangePoly[n][0], noNodes1D, F0);
8217 for(
MInt n = 0; n < nDim; n++) {
8219 &interp.
m_wBary[0], &lagrangePoly[n][0]);
8222 fill_n(state, noVars, 0.0);
8224 IF_CONSTEXPR(nDim == 2) {
8225 for(
MInt i = 0; i < noNodes1D; i++) {
8226 for(
MInt j = 0; j < noNodes1D; j++) {
8227 for(
MInt k = 0; k < noNodes1D3; k++) {
8228 for(
MInt v = 0; v < noVars; v++) {
8229 state[v] += U(i, j, k, v) * lagrangePoly[0][i] * lagrangePoly[1][j];
8236 for(
MInt i = 0; i < noNodes1D; i++) {
8237 for(
MInt j = 0; j < noNodes1D; j++) {
8238 for(
MInt k = 0; k < noNodes1D3; k++) {
8239 for(
MInt v = 0; v < noVars; v++) {
8240 state[v] += U(i, j, k, v) * lagrangePoly[0][i] * lagrangePoly[1][j] * lagrangePoly[2][k];
8278template <MInt nDim,
class SysEqn>
8280 const MInt elementId,
8282 calcStateAtPoint(point, elementId, SysEqn::noVars(), &m_elements.variables(elementId), state);
8287template <MInt nDim,
class SysEqn>
8290 const MBool NotUsed(interpolate)) {
8291 switch(sampleVarId) {
8294 calcStateAtPoint(point, elementId, state);
8299 calcStateAtPoint(point, elementId, SysEqn::noNodeVars(), &m_elements.nodeVars(elementId), state);
8304 calcStateAtPoint(point, elementId, SysEqn::noVars(), &m_elements.externalSource(elementId), state);
8308 TERMM(1,
"sampling variable not supported");
8331template <MInt nDim,
class SysEqn>
8334 RECORD_TIMER_START(m_timers[Timers::InitMortarProjection]);
8336 m_log <<
"Initializing mortar projections... ";
8338 using namespace dg::interpolation;
8345 m_projectionMatricesH.clear();
8346 m_projectionMatrixPointersH.clear();
8349 for(
MInt noNodes1D = m_minNoNodes1D; noNodes1D <= m_maxNoNodes1D; noNodes1D++) {
8350 sizeH += 4 * noNodes1D * noNodes1D;
8352 m_projectionMatricesH.resize(sizeH);
8355 m_projectionMatrixPointersH.resize(4 * (m_maxNoNodes1D + 1));
8357 for(
MInt noNodes1D = m_minNoNodes1D; noNodes1D <= m_maxNoNodes1D; noNodes1D++) {
8358 for(
MInt p = 0; p < 4; p++) {
8359 m_projectionMatrixPointersH[4 * noNodes1D + p] = &m_projectionMatricesH[offset];
8360 offset += noNodes1D * noNodes1D;
8366 for(
MInt noNodes1D = m_minNoNodes1D; noNodes1D <= m_maxNoNodes1D; noNodes1D++) {
8367 using namespace sbp::mortar;
8369 calcMortarProjectionMatrixHSBP(m_initPolyDeg, m_sbpOperator, mortarH<dg::mortar::forward>(noNodes1D, lower),
8370 mortarH<dg::mortar::forward>(noNodes1D, upper),
8371 mortarH<dg::mortar::reverse>(noNodes1D, lower),
8372 mortarH<dg::mortar::reverse>(noNodes1D, upper));
8375 for(
MInt polyDeg = m_minPolyDeg; polyDeg <= m_maxPolyDeg; polyDeg++) {
8376 const MInt noNodes1D = polyDeg + 1;
8380 using namespace dg::mortar;
8405 m_projectionMatricesP.clear();
8406 m_projectionMatrixPointersP.clear();
8409 for(
MInt polyDegHi = m_minPolyDeg; polyDegHi <= m_maxPolyDeg; polyDegHi++) {
8410 for(
MInt polyDegLo = m_minPolyDeg; polyDegLo < polyDegHi; polyDegLo++) {
8411 sizeP += 2 * (polyDegHi + 1) * (polyDegLo + 1);
8414 m_projectionMatricesP.resize(sizeP);
8417 m_projectionMatrixPointersP.resize(m_maxPolyDeg * (m_maxPolyDeg + 1));
8419 for(
MInt polyDegHi = m_minPolyDeg; polyDegHi <= m_maxPolyDeg; polyDegHi++) {
8420 for(
MInt polyDegLo = m_minPolyDeg; polyDegLo < polyDegHi; polyDegLo++) {
8421 for(
MInt p = 0; p < 2; p++) {
8422 const MInt idx = 2 * (polyDegHi * (polyDegHi - 1) / 2 + polyDegLo) + p;
8423 m_projectionMatrixPointersP[idx] = &m_projectionMatricesP[sizeP];
8424 sizeP += (polyDegHi + 1) * (polyDegLo + 1);
8432 for(
MInt polyDegHi = m_minPolyDeg; polyDegHi <= m_maxPolyDeg; polyDegHi++) {
8433 for(
MInt polyDegLo = m_minPolyDeg; polyDegLo < polyDegHi; polyDegLo++) {
8434 auto operatorLo = m_sbpOperator;
8435 auto operatorHi = m_sbpOperator;
8437 for(
MUint j = 0; j < m_prefPatchesPolyDeg.size(); j++) {
8438 if(polyDegLo == (
MInt)m_prefPatchesPolyDeg[j]) {
8439 operatorLo = m_prefPatchesOperators[j];
8441 if(polyDegHi == (
MInt)m_prefPatchesPolyDeg[j]) {
8442 operatorHi = m_prefPatchesOperators[j];
8447 mortarP<dg::mortar::forward>(polyDegLo, polyDegHi),
8448 mortarP<dg::mortar::reverse>(polyDegLo, polyDegHi));
8454 for(
MInt polyDegHi = m_minPolyDeg; polyDegHi <= m_maxPolyDeg; polyDegHi++) {
8455 const MInt noNodes1DHi = polyDegHi + 1;
8456 const DgInterpolation& interpHi = m_interpolation[polyDegHi][noNodes1DHi];
8457 for(
MInt polyDegLo = m_minPolyDeg; polyDegLo < polyDegHi; polyDegLo++) {
8458 const MInt noNodes1DLo = polyDegLo + 1;
8459 const DgInterpolation& interpLo = m_interpolation[polyDegLo][noNodes1DLo];
8469 mortarP<dg::mortar::forward>(polyDegLo, polyDegHi));
8474 mortarP<dg::mortar::reverse>(polyDegLo, polyDegHi));
8478 m_log <<
"done" << endl;
8480 RECORD_TIMER_STOP(m_timers[Timers::InitMortarProjection]);
8496template <MInt nDim,
class SysEqn>
8497template <MBool forward, MInt noVars>
8504 calcMortarProjectionH<forward, noVars>(srfcId, dir, source, destination, elem, surf);
8505 calcMortarProjectionP<forward, noVars>(srfcId, dir, source, destination, elem, surf);
8525template <MInt nDim,
class SysEqn>
8526template <MBool forward, MInt noVars>
8536 const MInt side = 1 - dir % 2;
8538 const MInt elementPolyDeg = elem.
polyDeg(elementId);
8540 if(elementPolyDeg == surfacePolyDeg) {
8546 const MInt surfaceNodes1D3 = (nDim == 3) ? surfaceNodes1D : 1;
8548 const MInt elementNodes1D3 = (nDim == 3) ? elementNodes1D : 1;
8554 src =
MFloatTensor(source, elementNodes1D, elementNodes1D3, noVars);
8555 dest =
MFloatTensor(destination, surfaceNodes1D, surfaceNodes1D3, noVars);
8557 src =
MFloatTensor(source, surfaceNodes1D, surfaceNodes1D3, noVars);
8558 dest =
MFloatTensor(destination, surfaceNodes1D, surfaceNodes1D3, noVars);
8565 const MInt size0 = forward ? surfaceNodes1D : elementNodes1D;
8566 const MInt size1 = forward ? elementNodes1D : surfaceNodes1D;
8569 MFloatMatrix p(mortarP<forward>(elementPolyDeg, surfacePolyDeg), size0, size1);
8573 IF_CONSTEXPR(nDim == 2) {
8575 for(
MInt i = 0; i < size0; i++) {
8576 for(
MInt j = 0; j < size1; j++) {
8577 for(
MInt var = 0; var < noVars; var++) {
8578 dest(i, 0, var) += src(j, 0, var) * p(i, j);
8585 for(
MInt i = 0; i < size0; i++) {
8586 for(
MInt s = 0; s < size0; s++) {
8587 for(
MInt j = 0; j < size1; j++) {
8588 for(
MInt k = 0; k < size1; k++) {
8589 for(
MInt var = 0; var < noVars; var++) {
8590 dest(i, s, var) += src(j, k, var) * p(i, j) * p(s, k);
8616template <MInt nDim,
class SysEqn>
8617template <MBool forward, MInt noVars>
8627 const MInt side = 1 - dir % 2;
8634 const MInt elementPolyDeg = elem.
polyDeg(elementId);
8638 const MInt noNodes1D = forward ? elementNoNodes1D : surfaceNoNodes1D;
8639 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
8642 MFloatTensor src(source, noNodes1D, noNodes1D3, noVars);
8643 MFloatTensor dest(destination, noNodes1D, noNodes1D3, noVars);
8654 const MInt dirA = (orientation == 0) ? 1 : 0;
8655 const MInt positionA = (grid().tree().coordinate(cellIdR, dirA) < grid().tree().coordinate(cellIdL, dirA))
8659 const MFloatMatrix pA(mortarH<forward>(noNodes1D, positionA), noNodes1D, noNodes1D);
8662 IF_CONSTEXPR(nDim == 2) {
8664 for(
MInt i = 0; i < noNodes1D; i++) {
8665 for(
MInt j = 0; j < noNodes1D; j++) {
8666 for(
MInt var = 0; var < noVars; var++) {
8667 dest(i, 0, var) += src(j, 0, var) * pA(i, j);
8675 const MInt dirB = (orientation == 2) ? 1 : 2;
8676 const MInt positionB = (grid().tree().coordinate(cellIdL, dirB) < grid().tree().coordinate(cellIdR, dirB))
8680 const MFloatMatrix pB(mortarH<forward>(noNodes1D, positionB), noNodes1D, noNodes1D);
8682 for(
MInt i = 0; i < noNodes1D; i++) {
8683 for(
MInt j = 0; j < noNodes1D; j++) {
8684 for(
MInt k = 0; k < noNodes1D; k++) {
8685 for(
MInt l = 0; l < noNodes1D; l++) {
8686 for(
MInt var = 0; var < noVars; var++) {
8687 dest(i, j, var) += src(k, l, var) * pA(i, k) * pB(j, l);
8697 if(elementPolyDeg != surfacePolyDeg) {
8699 copy_n(&dest[0], size, &src[0]);
8718template <MInt nDim,
class SysEqn>
8719template <MBool forward>
8722 ASSERT(sourcePolyDeg < targetPolyDeg,
8723 "source polynomial degree (= " + to_string(sourcePolyDeg)
8724 +
") must be lower than target polynomial degree (= " + to_string(targetPolyDeg) +
")");
8727 const MInt index = 2 * (targetPolyDeg * (targetPolyDeg - 1) / 2 + sourcePolyDeg) + (1 - forward);
8730 return m_projectionMatrixPointersP[index];
8747template <MInt nDim,
class SysEqn>
8748template <MBool forward>
8751 ASSERT(position == 0 || position == 1,
"position must be `0` or `1` (=" + to_string(position) +
")");
8754 const MInt index = 4 * noNodes1D + 2 * (1 - forward) + position;
8757 return m_projectionMatrixPointersH[index];
8766template <MInt nDim,
class SysEqn>
8770 m_log <<
"Exchanging polynomial degrees of MPI surfaces... ";
8773 vector<vector<MInt>> sendBuffer(m_noExchangeNghbrDomains);
8774 vector<vector<MInt>> recvBuffer(m_noExchangeNghbrDomains);
8775 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
8776 sendBuffer[i].resize(m_mpiSurfaces[i].size());
8777 recvBuffer[i].resize(m_mpiSurfaces[i].size());
8781 vector<MPI_Request> recvRequests(m_noExchangeNghbrDomains, MPI_REQUEST_NULL);
8782 vector<MPI_Request> sendRequests(m_noExchangeNghbrDomains, MPI_REQUEST_NULL);
8783 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
8785 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
8786 sendBuffer[i][j] = m_surfaces.polyDeg(m_mpiSurfaces[i][j]);
8791 m_exchangeNghbrDomains[i], mpiComm(), &recvRequests[i], AT_,
"recvBuffer[i][0]");
8793 domainId(), mpiComm(), &sendRequests[i], AT_,
"sendBuffer[i][0]");
8797 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
8798 MPI_Wait(&sendRequests[i], MPI_STATUS_IGNORE, AT_);
8799 MPI_Wait(&recvRequests[i], MPI_STATUS_IGNORE, AT_);
8803 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
8804 for(vector<MInt>::size_type j = 0; j < m_mpiSurfaces[i].size(); j++) {
8805 const MInt srfcId = m_mpiSurfaces[i][j];
8806 const MInt currentPolyDeg = m_surfaces.polyDeg(srfcId);
8807 const MInt receivedPolyDeg = recvBuffer[i][j];
8811 if(currentPolyDeg < receivedPolyDeg) {
8812 m_surfaces.polyDeg(srfcId) = receivedPolyDeg;
8813 m_surfaces.noNodes1D(srfcId) = receivedPolyDeg + 1;
8815 calcSurfaceNodeCoordinates(srfcId);
8820 m_log <<
"done" << endl;
8837template <MInt nDim,
class SysEqn>
8840 RECORD_TIMER_START(m_timers[Timers::AdaptiveRefinement]);
8843 set<MInt> adaptedSrfcs;
8846 const MInt noElements = m_elements.size();
8847 vector<MFloat> errorEstimate(noElements, 0.0);
8850 calcErrorEstimate(errorEstimate);
8853 MFloat localErrorEstimate = *max_element(errorEstimate.begin(), errorEstimate.end());
8857 if(noDomains() > 1) {
8859 "localErrorEstimate",
"buffer");
8862 const MFloat maxErrorEstimate = (buffer > localErrorEstimate) ? buffer : localErrorEstimate;
8865#pragma omp parallel for
8868 for(
MInt elementId = 0; elementId < noElements; elementId++) {
8869 const MInt cellId = m_elements.cellId(elementId);
8870 const MFloat error = errorEstimate[elementId];
8873 MInt adaptedPolyDeg = m_elements.polyDeg(elementId);
8880 if(grid().tree().coordinate(cellId, 0) > 0) {
8881 adaptedPolyDeg += 1;
8906 MFloat lvlThreshold = maxErrorEstimate * m_adaptiveThreshold;
8910 for(
MInt refLvl = m_maxPolyDeg; refLvl > m_minPolyDeg; refLvl--) {
8912 if(error > lvlThreshold || refLvl == m_minPolyDeg) {
8913 adaptedPolyDeg = refLvl;
8917 lvlThreshold *= m_adaptiveThreshold;
8920 TERMM(1,
"Invalid adaptive refinement case.");
8924 if(adaptedPolyDeg == m_elements.polyDeg(elementId)) {
8929 if(adaptedPolyDeg > m_maxPolyDeg) {
8930 m_log <<
"WARNING: Element polynomial degree would exceed maximum "
8931 "polynomial degree."
8933 cout <<
"WARNING: Element polynomial degree would exceed maximum "
8934 "polynomial degree."
8936 adaptedPolyDeg = m_maxPolyDeg;
8940 if(adaptedPolyDeg < m_minPolyDeg) {
8941 m_log <<
"WARNING: Element polynomial degree would be coarsed below "
8942 "minimum polynomial degree."
8944 cout <<
"WARNING: Element polynomial degree would be coarsed below "
8945 "minimum polynomial degree."
8947 adaptedPolyDeg = m_minPolyDeg;
8951 interpolateElement(elementId, adaptedPolyDeg);
8954 const MInt begin = 0;
8955 const MInt end = m_surfaces.size();
8958 MInt noAdaptedSrfcs = 0;
8960#pragma omp parallel for reduction(+ : noAdaptedSrfcs)
8962 for(
MInt srfcId = begin; srfcId < end; srfcId++) {
8963 const MInt srfcPolyDeg = m_surfaces.polyDeg(srfcId);
8964 const MInt internalSide = m_surfaces.internalSideId(srfcId);
8965 const MInt elementIdL = m_surfaces.nghbrElementIds(srfcId, (internalSide == -1) ? 0 : internalSide);
8966 const MInt elementIdR = m_surfaces.nghbrElementIds(srfcId, (internalSide == -1) ? 1 : internalSide);
8967 const MInt polyDegL = m_elements.polyDeg(elementIdL);
8968 const MInt polyDegR = m_elements.polyDeg(elementIdR);
8969 const MInt noNodes1DL = m_elements.noNodes1D(elementIdL);
8970 const MInt noNodes1DR = m_elements.noNodes1D(elementIdR);
8974 if(max(polyDegL, polyDegR) != srfcPolyDeg) {
8975 m_surfaces.polyDeg(srfcId) = max(polyDegL, polyDegR);
8976 m_surfaces.noNodes1D(srfcId) = max(noNodes1DL, noNodes1DR);
8977 calcSurfaceNodeCoordinates(srfcId);
8983 if(hasMpiExchange()) {
8984 exchangeMpiSurfacePolyDeg();
8988 m_noTotalNodesXD = 0;
8989 for(
MInt i = 0; i < m_elements.size(); i++) {
8990 const MInt noNodesXD = m_elements.noNodesXD(i);
8991 m_noTotalNodesXD += noNodesXD;
8994 cout <<
"Grid has been refined at timestep " << timeStep <<
" (adapted surfaces: " << noAdaptedSrfcs <<
" )" << endl;
8995 m_log <<
"Grid has been refined at timestep " << timeStep <<
" (adapted surfaces: " << noAdaptedSrfcs <<
" )" << endl;
8997 RECORD_TIMER_STOP(m_timers[Timers::AdaptiveRefinement]);
9008template <MInt nDim,
class SysEqn>
9015 const MInt begin = m_innerSurfacesOffset;
9016 const MInt end = m_surfaces.size();
9017 const MInt noVars = SysEqn::noVars();
9018 const MInt noElements = m_elements.size();
9021 fill_n(&error[0], m_surfaces.
size(), 0.0);
9025#pragma omp parallel for
9027 for(
MInt srfcId = begin; srfcId < end; srfcId++) {
9028 const MInt noNodes1D = m_surfaces.noNodes1D(srfcId);
9029 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
9030 const MFloat* stateL = &m_surfaces.variables(srfcId, 0);
9031 const MFloat* stateR = &m_surfaces.variables(srfcId, 1);
9038 for(
MInt i = 0; i < noNodes1D; i++) {
9039 for(
MInt j = 0; j < noNodes1D3; j++) {
9040 for(
MInt var = 0; var < noVars; var++) {
9041 error[srfcId] += fabs(uL(i, j, var) - uR(i, j, var));
9047 error[srfcId] = error[srfcId] / (noNodes1D * noNodes1D3);
9051#pragma omp parallel for
9054 for(
MInt elementId = 0; elementId < noElements; elementId++) {
9055 const MInt noSrfcs = 2 * nDim;
9056 for(
MInt i = 0; i < noSrfcs; i++) {
9057 const MInt srfcId = m_elements.surfaceIds(elementId, i);
9058 errorEstimate[elementId] += error[srfcId];
9074template <MInt nDim,
class SysEqn>
9081 const MInt adaptedNodes1D = adaptedPolyDeg + 1;
9082 const MInt adaptedNodes1D3 = (nDim == 3) ? adaptedNodes1D : 1;
9083 const MInt noVars = SysEqn::noVars();
9084 const MInt elemPolyDeg = m_elements.polyDeg(elementId);
9085 const MInt elemNodes1D = m_elements.noNodes1D(elementId);
9086 const MInt size = adaptedNodes1D * adaptedNodes1D * adaptedNodes1D3 * noVars;
9087 vector<MFloat> buffer(size);
9088 using namespace dg::interpolation;
9091 if(elemPolyDeg < adaptedPolyDeg) {
9092 interpolateNodes<nDim>(&m_elements.variables(elementId), mortarP<dg::mortar::forward>(elemPolyDeg, adaptedPolyDeg),
9093 elemNodes1D, adaptedNodes1D, noVars, &buffer[0]);
9094 copy(buffer.begin(), buffer.end(), &m_elements.variables(elementId));
9097 interpolateNodes<nDim>(&m_elements.timeIntStorage(elementId),
9098 mortarP<dg::mortar::forward>(elemPolyDeg, adaptedPolyDeg), elemNodes1D, adaptedNodes1D,
9099 noVars, &buffer[0]);
9100 copy(buffer.begin(), buffer.end(), &m_elements.timeIntStorage(elementId));
9102 interpolateNodes<nDim>(&m_elements.variables(elementId), mortarP<dg::mortar::reverse>(adaptedPolyDeg, elemPolyDeg),
9103 elemNodes1D, adaptedNodes1D, noVars, &buffer[0]);
9104 copy(buffer.begin(), buffer.end(), &m_elements.variables(elementId));
9107 interpolateNodes<nDim>(&m_elements.timeIntStorage(elementId),
9108 mortarP<dg::mortar::reverse>(adaptedPolyDeg, elemPolyDeg), elemNodes1D, adaptedNodes1D,
9109 noVars, &buffer[0]);
9110 copy(buffer.begin(), buffer.end(), &m_elements.timeIntStorage(elementId));
9114 m_elements.polyDeg(elementId) = adaptedPolyDeg;
9115 m_elements.noNodes1D(elementId) = adaptedNodes1D;
9118 calcElementNodeCoordinates(elementId);
9135template <MInt nDim,
class SysEqn>
9139 const MInt cellId = m_elements.cellId(elementId);
9140 MInt nghbrCellId = grid().tree().neighbor(cellId, dir);
9143 if(nghbrCellId < 0 && grid().tree().hasParent(cellId)) {
9144 nghbrCellId = grid().tree().neighbor(grid().tree().parent(cellId), dir);
9147 return nghbrCellId >= 0;
9158template <MInt nDim,
class SysEqn>
9161 ASSERT(elementId >= 0,
"Invalid elementId");
9163 const MInt cellId = m_elements.cellId(elementId);
9164 const MInt level = grid().tree().level(cellId);
9204template <MInt nDim,
class SysEqn>
9208 const MInt nghbrCellId = grid().tree().neighbor(cellId, dir);
9209 const MInt cellLvl = grid().tree().level(cellId);
9212 if(nghbrCellId > 0 && grid().tree().hasChildren(nghbrCellId)) {
9213 TERMM(1,
"The following is not yet tested! If it works, delete this TERMM!");
9214 for(
MInt child = 0; child <
IPOW2(nDim); child++) {
9215 if(!(childCodePro[dir] & (1 << child)))
continue;
9216 if(grid().tree().child(nghbrCellId, child) > -1)
return cellLvl + 1;
9221 if(nghbrCellId < 0 && grid().tree().hasParent(cellId)) {
9222 if(grid().tree().neighbor(grid().tree().parent(cellId), dir)) {
9227 return nghbrCellId < 0 ? -1 : cellLvl;
9234template <MInt nDim,
class SysEqn>
9238 return (m_adaptivePref > 0 && hasPref());
9244template <MInt nDim,
class SysEqn>
9248 return (m_pref == 1 && m_minPolyDeg != m_maxPolyDeg);
9258template <MInt nDim,
class SysEqn>
9265 return (hasAdaptivePref() && (timeStep == 2 || timeStep % m_adaptiveInterval == 0));
9275template <MInt nDim,
class SysEqn>
9279 const MInt begin = m_mpiSurfacesOffset;
9280 const MInt end = m_mpiSurfacesOffset + m_noMpiSurfaces;
9281 MBool isMpiSrfc =
false;
9283 if(begin <=
id &&
id < end) {
9297template <MInt nDim,
class SysEqn>
9301 MInt hElementId = -1;
9303 const MInt noHElements = m_helements.size();
9304 for(
MInt hId = 0; hId < noHElements; hId++) {
9305 if(m_helements.elementId(hId) == elementId) {
9318template <MInt nDim,
class SysEqn>
9322 finalizeMpiExchange();
9325 m_maxNoSurfaces = -1;
9328 for(
MInt i = 0; i < m_noExchangeNghbrDomains; i++) {
9329 m_mpiSurfaces[i].clear();
9330 m_sendBuffers[i].clear();
9331 m_recvBuffers[i].clear();
9333 m_mpiSurfaces.clear();
9334 m_sendBuffers.clear();
9335 m_recvBuffers.clear();
9336 m_exchangeNghbrDomains.clear();
9337 m_sendRequests.clear();
9338 m_recvRequests.clear();
9339 m_noExchangeNghbrDomains = 0;
9341 m_isInitSolver =
false;
9348template <MInt nDim,
class SysEqn>
9352 MInt totalNoDgLT = 0;
9355 totalNoDgLT += m_maxPolyDeg - m_minPolyDeg + 1;
9358 if(m_weightDgRbcElements) {
9367template <MInt nDim,
class SysEqn>
9371 const MInt noPolyDegs = m_maxPolyDeg - m_minPolyDeg + 1;
9372 for(
MInt i = 0; i < noPolyDegs; i++) {
9374 names[i] =
"dg_node_p" + std::to_string(m_minPolyDeg + i);
9377 if(m_weightDgRbcElements) {
9378 weights[noPolyDegs] = 0.5;
9379 names[noPolyDegs] =
"dg_node_rbc";
9389template <MInt nDim,
class SysEqn>
9394 std::fill_n(&loadQuantities[0], noLoadTypes(), 0);
9401 const MInt noPolyDegs = m_maxPolyDeg - m_minPolyDeg + 1;
9403 if(noPolyDegs > 1) {
9405 for(
MInt i = 0; i < noPolyDegs; i++) {
9406 loadQuantities[i] = m_statLocalNoActiveDOFsPolyDeg[i];
9410 loadQuantities[0] = m_statLocalNoActiveDOFs;
9414 if(m_weightDgRbcElements) {
9415 loadQuantities[noPolyDegs] = 0;
9416 for(
auto&& bc : m_boundaryConditions) {
9417 loadQuantities[noPolyDegs] += bc->getLocalNoNodes();
9430template <MInt nDim,
class SysEqn>
9433 ASSERT(isActive(),
"solver is not active");
9436 const MInt cellId = grid().tree().grid2solver(gridCellId);
9441 if(cellId < 0 || cellId >= grid().noInternalCells()) {
9442 TERMM(1,
"The given cell id is invalid.");
9445 const MInt elementId = m_elements.getElementByCellId(cellId);
9451 if(elementId != -1) {
9452 const MInt polyDeg = m_elements.polyDeg(elementId);
9453 const MInt noNodesXD = m_elements.noNodesXD(elementId);
9456 cellLoad = m_weightPerElement + noNodesXD * weights[polyDeg - m_minPolyDeg];
9459 if(m_weightDgRbcElements) {
9460 for(
auto&& bc : m_boundaryConditions) {
9461 if(bc->hasBcElement(elementId)) {
9462 cellLoad += noNodesXD * weights[m_maxPolyDeg - m_minPolyDeg + 1];
9473template <MInt nDim,
class SysEqn>
9476 ASSERT(isActive(),
"solver is not active");
9478 const MInt noCells = grid().noInternalCells();
9479 const MInt noCellsGrid = grid().raw().treeb().size();
9480 const MInt noWeights = noLoadTypes();
9481 const MInt offset = solverId() * noCellsGrid;
9483 std::vector<MFloat> weights(noWeights);
9484 std::fill(weights.begin(), weights.end(), m_weightPerNode);
9486 for(
MInt cellId = 0; cellId < noCells; cellId++) {
9487 const MInt gridCellId = grid().tree().solver2grid(cellId);
9488 solverCellWeight[offset + gridCellId] = getCellLoad(gridCellId, &weights[0]);
9496template <MInt nDim,
class SysEqn>
9505 const MInt domainOffset = grid().domainOffset(domainId());
9506 const MInt noElements = m_elements.size();
9508 for(
MInt elementId = 0; elementId < noElements; elementId++) {
9509 const MInt globalCellId = m_elements.cellId(elementId) + domainOffset;
9510 m_elements.cellId(elementId) = globalCellId;
9518template <MInt nDim,
class SysEqn>
9527 const MInt domainOffset = grid().domainOffset(domainId());
9528 const MInt noElements = m_elements.size();
9530 for(
MInt elementId = 0; elementId < noElements; elementId++) {
9531 const MInt localCellId = m_elements.cellId(elementId) - domainOffset;
9532 m_elements.cellId(elementId) = localCellId;
9543template <MInt nDim,
class SysEqn>
9548 TERMM(1,
"Error: cellDataTypeDlb() might give wrong results on inactive ranks.");
9553 if(dataId > -1 && dataId < noDgCartesianSolverCellData()) {
9555 dataType = s_cellDataTypeDlb[dataId];
9556 }
else if(dataId >= noDgCartesianSolverCellData() && dataId < noCellDataDlb()) {
9558 MInt offset = noDgCartesianSolverCellData();
9559 for(
auto&& bc : m_boundaryConditions) {
9560 const MInt bcNoCellData = bc->noCellDataDlb();
9561 if(dataId >= offset && dataId < offset + bcNoCellData) {
9562 dataType = bc->cellDataTypeDlb(dataId - offset);
9565 offset += bcNoCellData;
9568 TERMM(1,
"The requested dataId is not valid: " + to_string(dataId) +
" (" + to_string(noDgCartesianSolverCellData())
9569 +
", " + to_string(noCellDataDlb()) +
")");
9583template <MInt nDim,
class SysEqn>
9593 const MInt cellId = grid().tree().grid2solver(gridCellId);
9600 if(dataId > -1 && dataId < noDgCartesianSolverCellData()) {
9602 const MInt elementId = m_elements.getElementByCellId(cellId);
9604 if(elementId != -1) {
9605 const MInt noNodesXD = m_elements.noNodesXD(elementId);
9607 case CellData::ELEM_CELL_ID:
9608 case CellData::ELEM_POLY_DEG:
9611 case CellData::ELEM_NO_NODES_1D:
9614 case CellData::ELEM_VARIABLES:
9615 dataSize = noNodesXD * SysEqn::noVars();
9617 case CellData::ELEM_NODE_VARS:
9618 dataSize = noNodesXD * SysEqn::noNodeVars();
9621 TERMM(1,
"Unknown data id. (" + to_string(dataId) +
")");
9625 }
else if(dataId >= noDgCartesianSolverCellData() && dataId < noCellDataDlb()) {
9627 MInt offset = noDgCartesianSolverCellData();
9628 for(
auto&& bc : m_boundaryConditions) {
9629 const MInt bcNoCellData = bc->noCellDataDlb();
9630 if(dataId >= offset && dataId < offset + bcNoCellData) {
9631 dataSize = bc->cellDataSizeDlb(dataId - offset, cellId);
9634 offset += bcNoCellData;
9637 TERMM(1,
"The requested dataId is not valid.");
9647template <MInt nDim,
class SysEqn>
9652 m_loadBalancingReinitStage = 0;
9665 updateDomainInfo(-1, -1, MPI_COMM_NULL, AT_);
9666 m_elements.reset(0);
9667 m_helements.reset(0);
9668 m_surfaces.reset(0);
9673 updateDomainInfo(grid().domainId(), grid().noDomains(), grid().mpiComm(), AT_);
9675 allocateAndInitSolverMemory();
9678 printAllocatedMemory(previouslyAllocated,
"DgCartesianSolver (solverId = " + to_string(m_solverId) +
")", mpiComm());
9680 RECORD_TIMER_START(m_timers[Timers::RunInit]);
9681 RECORD_TIMER_START(m_timers[Timers::InitSolverObjects]);
9687 RECORD_TIMER_STOP(m_timers[Timers::InitSolverObjects]);
9688 RECORD_TIMER_STOP(m_timers[Timers::RunInit]);
9695template <MInt nDim,
class SysEqn>
9705 m_loadBalancingReinitStage = 1;
9707 RECORD_TIMER_START(m_timers[Timers::RunInit]);
9708 RECORD_TIMER_START(m_timers[Timers::InitSolverObjects]);
9710 initInterpolation();
9713 m_internalDataSize = m_elements.size() * pow(m_maxNoNodes1D, nDim) * SysEqn::noVars();
9716 m_noTotalNodesXD = 0;
9717 for(
MInt i = 0; i < m_elements.size(); i++) {
9718 const MInt noNodesXD = m_elements.noNodesXD(i);
9719 m_noTotalNodesXD += noNodesXD;
9722 initNodeCoordinates();
9726 checkCellProperties();
9733 if(noDomains() > 1) {
9738 m_log <<
"Reinitializing sponge... ";
9739 sponge().init(m_maxPolyDeg, &grid(), &m_elements, &m_surfaces, &m_boundaryConditions, &m_sysEqn, mpiComm());
9740 m_log <<
"done" << endl;
9743 if(noDomains() > 1 && hasPref()) {
9744 exchangeMpiSurfacePolyDeg();
9747 initMortarProjections();
9752 m_surfaceData.init();
9753 m_volumeData.init();
9757 if(m_pointData.enabled() || m_surfaceData.enabled()) {
9758 TERMM(1,
"DLB for point/surface/volumeData not supported yet");
9761 initSimulationStatistics();
9768 if(SysEqn::noNodeVars() > 0) {
9769 updateNodeVariables();
9770 extendNodeVariables();
9773 m_isInitSolver =
true;
9774 m_isInitData =
true;
9776 m_isInitMainLoop =
true;
9779 RECORD_TIMER_STOP(m_timers[Timers::InitSolverObjects]);
9780 RECORD_TIMER_STOP(m_timers[Timers::RunInit]);
9783 m_loadBalancingReinitStage = 2;
9788template <MInt nDim,
class SysEqn>
9790 std::vector<MInt>& globalIntVars) {
9793 globalFloatVars.push_back(m_time);
9794 globalFloatVars.push_back(m_dt);
9796 globalIntVars.push_back(m_timeStep);
9797 globalIntVars.push_back(m_firstTimeStep);
9798 globalIntVars.push_back(m_noAnalyzeTimeSteps);
9803template <MInt nDim,
class SysEqn>
9805 std::vector<MInt>& globalIntVars) {
9808 m_time = globalFloatVars[0];
9809 m_dt = globalFloatVars[1];
9811 m_timeStep = globalIntVars[0];
9812 m_firstTimeStep = globalIntVars[1];
9813 m_noAnalyzeTimeSteps = globalIntVars[2];
9818template <MInt nDim,
class SysEqn>
9820 const MBool allTimings) {
9822 const MString namePrefix =
"b" + std::to_string(solverId()) +
"_";
9824 const MFloat load = returnLoadRecord();
9825 const MFloat idle = returnIdleRecord();
9827 solverTimings.emplace_back(namePrefix +
"loadDgCartesianSolver", load);
9828 solverTimings.emplace_back(namePrefix +
"idleDgCartesianSolver", idle);
9830#ifdef MAIA_TIMER_FUNCTION
9831 solverTimings.emplace_back(namePrefix +
"timeStepRk", RETURN_TIMER_TIME(m_timers[Timers::RungeKuttaStep]));
9835 solverTimings.emplace_back(namePrefix +
"calcDgTimeDerivative", RETURN_TIMER_TIME(m_timers[Timers::TimeDeriv]));
9837 solverTimings.emplace_back(namePrefix +
"resetRHS", RETURN_TIMER_TIME(m_timers[Timers::ResetRHS]));
9838 solverTimings.emplace_back(namePrefix +
"prolong_to_surfaces", RETURN_TIMER_TIME(m_timers[Timers::Prolong]));
9839 solverTimings.emplace_back(namePrefix +
"forward_projection",
9840 RETURN_TIMER_TIME(m_timers[Timers::ForwardProjection]));
9842 solverTimings.emplace_back(namePrefix +
"MPI_surface_exchange", RETURN_TIMER_TIME(m_timers[Timers::SurfExchange]));
9843 solverTimings.emplace_back(namePrefix +
"communication", RETURN_TIMER_TIME(m_timers[Timers::SurfExchangeComm]));
9844 solverTimings.emplace_back(namePrefix +
"copy_operations", RETURN_TIMER_TIME(m_timers[Timers::SurfExchangeCopy]));
9845 solverTimings.emplace_back(namePrefix +
"waiting", RETURN_TIMER_TIME(m_timers[Timers::SurfExchangeWait]));
9846 solverTimings.emplace_back(namePrefix +
"waiting_send", RETURN_TIMER_TIME(m_timers[Timers::SEWaitSend]));
9847 solverTimings.emplace_back(namePrefix +
"waiting_recv", RETURN_TIMER_TIME(m_timers[Timers::SEWaitRecv]));
9849 solverTimings.emplace_back(namePrefix +
"calcVolumeIntegral", RETURN_TIMER_TIME(m_timers[Timers::VolInt]));
9851 solverTimings.emplace_back(namePrefix +
"flux_calculation", RETURN_TIMER_TIME(m_timers[Timers::Flux]));
9852 solverTimings.emplace_back(namePrefix +
"calcBoundarySurfaceFlux", RETURN_TIMER_TIME(m_timers[Timers::FluxBndry]));
9853 solverTimings.emplace_back(namePrefix +
"calcInnerSurfaceFlux", RETURN_TIMER_TIME(m_timers[Timers::FluxInner]));
9854 solverTimings.emplace_back(namePrefix +
"calcMpiSurfaceFlux", RETURN_TIMER_TIME(m_timers[Timers::FluxMPI]));
9856 solverTimings.emplace_back(namePrefix +
"surface_integrals", RETURN_TIMER_TIME(m_timers[Timers::SurfInt]));
9857 solverTimings.emplace_back(namePrefix +
"applyJacobian", RETURN_TIMER_TIME(m_timers[Timers::Jacobian]));
9858 solverTimings.emplace_back(namePrefix +
"calcSourceTerms", RETURN_TIMER_TIME(m_timers[Timers::Sources]));
9859 solverTimings.emplace_back(namePrefix +
"resetExtSources",
9860 RETURN_TIMER_TIME(m_timers[Timers::ResetExternalSources]));
9861 solverTimings.emplace_back(namePrefix +
"applyExternalSources",
9862 RETURN_TIMER_TIME(m_timers[Timers::ExternalSources]));
9863 solverTimings.emplace_back(namePrefix +
"calcSpongeTerms", RETURN_TIMER_TIME(m_timers[Timers::Sponge]));
9865 solverTimings.emplace_back(namePrefix +
"time_integration", RETURN_TIMER_TIME(m_timers[Timers::TimeInt]));
9867 solverTimings.emplace_back(namePrefix +
"IO", RETURN_TIMER_TIME(m_timers[Timers::MainLoopIO]));
9868 solverTimings.emplace_back(namePrefix +
"solution_analysis", RETURN_TIMER_TIME(m_timers[Timers::Analysis]));
9869 solverTimings.emplace_back(namePrefix +
"adaptive_refinement",
9870 RETURN_TIMER_TIME(m_timers[Timers::AdaptiveRefinement]));
9874 solverTimings.emplace_back(namePrefix +
"MPI_surface_exchange", RETURN_TIMER_TIME(m_timers[Timers::SurfExchange]));
9876 solverTimings.emplace_back(namePrefix +
"calcBoundarySurfaceFlux", RETURN_TIMER_TIME(m_timers[Timers::FluxBndry]));
9883template <MInt nDim,
class SysEqn>
9885 std::vector<std::pair<MString, MInt>>& domainInfo) {
9888 const MString namePrefix =
"b" + std::to_string(solverId()) +
"_";
9891 const MInt noElements = m_elements.size();
9892 domainInfo.emplace_back(namePrefix +
"noDgElements", noElements);
9895 MInt noBcElements = 0;
9896 for(
auto&& bc : m_boundaryConditions) {
9897 for(
MInt elementId = 0; elementId < noElements; elementId++) {
9898 if(bc->hasBcElement(elementId)) {
9903 domainInfo.emplace_back(namePrefix +
"noDgBcElements", noBcElements);
9908template <MInt nDim,
class SysEqn>
9910 std::vector<MInt>& noSamplingVars,
9911 std::vector<std::vector<MString>>& samplingVarNames,
9916 std::vector<MString> varNamesList;
9917 MInt noSampleVars = readSolverSamplingVarNames(varNamesList, featureName);
9920 if(noSampleVars == 0) {
9921 varNamesList.push_back(
"DG_VARS");
9925 for(
MInt i = 0; i < noSampleVars; i++) {
9927 std::vector<MString> varNames;
9929 auto samplingVarIt = std::find(samplingVarIds.begin(), samplingVarIds.end(), samplingVar);
9930 if(samplingVarIt != samplingVarIds.end()) {
9931 TERMM(1,
"Sampling variable '" + varNamesList[i] +
"' already specified.");
9934 switch(samplingVar) {
9936 const MInt noVars = SysEqn::noVars();
9938 samplingVarIds.push_back(
DG_VARS);
9939 noSamplingVars.push_back(noVars);
9941 varNames.resize(noVars);
9942 for(
MInt varId = 0; varId < noVars; varId++) {
9944 varNames[varId] = SysEqn::consVarNames(varId);
9947 samplingVarNames.push_back(varNames);
9951 const MInt noVars = SysEqn::noNodeVars();
9954 noSamplingVars.push_back(noVars);
9956 varNames.resize(noVars);
9957 for(
MInt varId = 0; varId < noVars; varId++) {
9958 varNames[varId] = SysEqn::nodeVarNames(varId);
9960 samplingVarNames.push_back(varNames);
9964 const MInt noVars = SysEqn::noVars();
9967 noSamplingVars.push_back(noVars);
9969 varNames.resize(noVars);
9970 for(
MInt varId = 0; varId < noVars; varId++) {
9971 varNames[varId] =
"source_" + SysEqn::consVarNames(varId);
9973 samplingVarNames.push_back(varNames);
9977 TERMM(1,
"Unknown sampling variable: " + varNamesList[i]);
9986#if defined(MAIA_GCC_COMPILER)
9987#pragma GCC diagnostic pop
MLong allocatedBytes()
Return the number of allocated bytes.
static MInt propertyLength(const MString &name, MInt solverId=m_noSolvers)
Returns the number of elements of a property.
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
MBool needHElementForCell(const MInt cellId)
Return true if h-element is needed for cell, false otherwise.
std::array< MInt, 2 *nDim > getBoundaryConditionIds(const MInt cellId)
Determine if element is boundary is cut by geometry elements and return corresponding boundary condit...
MInt noLoadTypes() const override
Return the number of DG load types.
void resetBuffer(const MInt totalSize, MFloat *const buffer)
Reset the given buffer to zero.
void resetExternalSources()
DgBoundaryConditionFactory< nDim, SysEqn > m_boundaryConditionFactory
MBool hasMpiExchange() const
void prolongToSurfaces()
Extrapolate the solution from inside the elements to the surfaces.
MFloat * mortarP(const MInt sourcePolyDeg, const MInt targetPolyDeg)
void checkCellProperties()
Check all relevant bit properties in the cells.
void calcSurfaceIntegral()
Calculate the surface integral for all faces of element and update dU/dt.
void applyJacobian()
Adds the negative of the inverse Jacobian to the time derivative.
void preTimeStep() override
Perform pre-time-step operations, e.g. set new dt if required.
void initTimers()
Initialize all solver-wide timers and start the solver timer.
void calcInnerSurfaceFlux()
Calculate the numerical flux on the internal surfaces and update m_flux.
MBool calcStateAtPoint(const MFloat *point, MFloat *state)
returns the cellId of a cell containing a given point
typename DgBoundaryConditionFactory< nDim, SysEqn >::ReturnType BC
void calcMortarProjectionH(const MInt srfcId, const MInt dir, MFloat *source, MFloat *destination, ElementCollector &elem, SurfaceCollector &surf)
Calculate the h-refinement mortar projection.
MInt getLevelOfDirectNeighborCell(const MInt cellId, const MInt dir)
void applyExternalSourceTerms(const MFloat time)
Add the external coupling source terms to the right hand side.
MInt getHElementId(const MInt elementId) const
Return h-element id of a given element (if it exists).
void createHElement(const MInt cellId)
Create h-element for cell with id cellId.
MFloat calcTimeStep()
Calculate the next time step.
void loadNodeVarsFile()
Load node variables from file.
void calcSamplingVarAtPoint(const MFloat *point, const MInt elementId, const MInt sampleVarId, MFloat *state, const MBool interpolate) override
Calculate the state vector at a given point and for the specified sampling variable.
void createElement(const MInt cellId)
Create element for cell with id cellId.
void initSurfaces()
Create for all elements and directions surfaces if necessary.
void setGlobalSolverVars(std::vector< MFloat > &NotUsed(globalFloatVars), std::vector< MInt > &NotUsed(globalIdVars)) override
Set global solver variables for DLB (same on each domain)
MBool pointIsInside(const MFloat *const coordinates)
Check if point is inside the computational domain.
void calcMpiSurfaceFlux()
Calculate the numerical flux on the MPI surfaces and update m_flux.
void initJacobian()
Calculates the inverse Jacobian for each element.
void exchangeMpiSurfacePolyDeg()
MBool hasAdaptivePref() const
Return true if adaptive hp-refinement is activated.
void getSolverSamplingProperties(std::vector< MInt > &samplingVarIds, std::vector< MInt > &noSamplingVars, std::vector< std::vector< MString > > &samplingVarNames, const MString featureName="") override
Return sampling properties for the DG solver.
void getLoadQuantities(MInt *const loadQuantities) const override
Return the cumulative load quantities on this domain.
void writeRestartFile(const MBool writeRestart, const MBool writeBackup, const MString gridFileName, MInt *recalcIdTree) override
Write a restart file.
MInt initMpiSurface(const MInt elementId, const MInt dir)
Check if the surface of the element in the given direction is a MPI surface. If it is,...
MInt calculateNeededNoSurfaces()
Determine the number of surfaces that need to be created on this domain.
MBool solutionStep() override
Perform one Runge-Kutta step/stage.
MInt initInnerSurface(const MInt elementId, const MInt dir)
Check if the surface of the element in the given direction is an inner surface without h/p-refinement...
void cleanUp() override
Clean up after a simulation run.
MString getRestartFileName(const MInt timeStep, const MInt useNonSpecifiedRestartFile)
Return name of a restart file for the given time step.
virtual void saveRestartFile()
Saves a file to disk with all information that is necessary to restart the calculations from here.
MBool isMpiSurface(const MInt id) const
Return true if a surface is a MPI surface.
MBool hasSurface(const MInt elementId, const MInt dir)
Check if a surface exists for the element in the given direction.
void calcMortarProjection(const MInt srfcId, const MInt dir, MFloat *source, MFloat *destination, ElementCollector &elem, SurfaceCollector &surf)
void initDgCartesianSolver()
Initializes basic properties of the solver.
void globalToLocalIds() override
Change global into local ids.
void initMortarProjections()
Calculate all necessary mortar projection matrices.
MInt createHMPISurfaces(const MInt elementId, const MInt dir)
void setInputOutputProperties()
MInt cellDataTypeDlb(const MInt dataId) const override
Get data type of cell data.
void checkGridMap(const MString &donorGridFileName)
Perform some sanity checks on loaded grid map.
void setCellWeights(MFloat *solverCellWeight) override
Set cell weights with constant weighting factor weightPerNode.
std::array< MInt, Timers::_count > m_timers
void calcSourceTerms(MFloat t)
Calculates the source terms for each node and adds them to the time derivative of the conservative va...
void saveNodeVarsFile()
Save node variables to file.
void extendNodeVariables()
Extends nodeVars from given planes to given directions.
void updateNodeVariablesSingleElement(const MInt elementId, const MInt extendDir, MIntTensor hForwardSurfaceId, MIntTensor hReverseSurfaceId, std::vector< MBool > &cellHasUpdatedMeans, std::vector< MInt > &SurfaceWantsToMPISend, MBool &noMoreUpdates)
Sets nodeVars of an element to values on the surface opposite to extendDir.
void adaptiveRefinement(const MInt timeStep)
Apply adaptive refinement (right now: only p-refinement)
void initNodeCoordinates()
Calculates the coordinates of the integration nodes within each element.
~DgCartesianSolver() override
Frees resources that were allocated in the constructor and that are not automatically released when t...
void calcVolumeIntegral()
Calculate the volume integral for all elements and update m_rightHandSide.
void calcElementNodeCoordinates(const MInt elementId)
Calculates the coordinates of the integration nodes within the element.
DgCartesianSolver(const MInt solverId, GridProxy &gridProxy_, Geometry< nDim > &geometry_, const MPI_Comm comm)
Constructor of the DG solver reads properties and allocates solver resources (as far as they are know...
MFloat * mortarH(const MInt polyDeg, const MInt position)
void initPolynomialDegree()
Calculate and set initial polynomial degree and number of nodes in all elements.
void loadGridMap(const MString &gridMapFileName)
Load previously created grid map.
void initInterpolation()
Calculates necessary coefficients for interpolation and stores them once for the whole solver.
void initElements()
Initialize all elements by iterating over all cells and creating an element for each internal leaf ce...
void calcMortarProjectionP(const MInt srfcId, const MInt dir, MFloat *source, MFloat *destination, ElementCollector &elem, SurfaceCollector &surf)
Calculate the p-refinement mortar projection.
void calcBoundarySurfaceFlux(MFloat t)
Calculate the numerical flux on the boundary surfaces and update m_flux.
void saveSolverSolution(const MBool forceOutput, const MBool finalTimeStep) override
Save the solver solution, i.e. write solution files and sampling/slice output.
void calcErrorEstimate(std::vector< MFloat > &errorEstimate)
Calculate error estimate for adaptive hp-refinement.
void initSolver() override
Initialize the solver.
void allocateAndInitSolverMemory()
Allocates main memory resources.
void analyzeTimeStep(MFloat t, MFloat runTimeRelative, MFloat runTimeTotal, MInt timeStep, MFloat dt)
Calculates and prints the L^2 and L^inf errors (using calcErrorNorms()) for the current time.
void initData()
Initialize solver data, i.e. set initial conditions or load restart file.
MFloat getCellLoad(const MInt cellId, const MFloat *const weights) const override
Return the load of a single cell (given computational weights).
MInt getLevelByElementId(const MInt elementId)
MBool needElementForCell(const MInt cellId)
Return true if element is needed for cell, false otherwise.
MBool hasPref() const
Return true if p-refinement is set.
MBool hasNeighborCell(const MInt elementId, const MInt dir)
void getSolverTimings(std::vector< std::pair< MString, MFloat > > &solverTimings, const MBool allTimings) override
Get solver timings.
MInt cellDataSizeDlb(const MInt dataId, const MInt cellId) override
Return data size of cell data.
void interpolateElement(const MInt elementId, const MInt newPolyDeg)
Interpolate an element to a different polynomial degree.
void initMainLoop()
Perform all operations that prepare the execution of the main loop.
void initSimulationStatistics()
void getGlobalSolverVars(std::vector< MFloat > &NotUsed(globalFloatVars), std::vector< MInt > &NotUsed(globalIntVars)) override
Get global solver variables that need to be communicated for DLB (same on each domain)
void getDomainDecompositionInformation(std::vector< std::pair< MString, MInt > > &domainInfo) override
Return decomposition information, i.e. number of local elements,...
void balancePost() override
Reinitialize solver after setting solution data.
void calcErrorNorms(const MFloat t, std::vector< MFloat > &L2Error, std::vector< MFloat > &LInfError, std::vector< MFloat > &L2ErrLocal, std::vector< MFloat > &LInfErrLocal)
Calculate the L^2 and L^infinity error norms at a given time. The error is calculated in comparison t...
typename CartesianSolver::GridProxy GridProxy
void applySurfaceIntegral(MFloat *rhs, const MInt polyDeg, const MInt noNodes1D, const MInt srfcId, const MInt side, const MFloat *flux, SurfaceCollector &surf)
Calculate the surface integral for a face of an element and update dU/dt.
MInt initBoundarySurface(const MInt elementId, const MInt dir)
Check if the surface of the element in the given direction is a boundary surface. If it is init,...
MInt createSurface(const MInt elementId, const MInt dir, MInt nghbrId=-1)
void cancelMpiRequests() override
Cancel open MPI (receive) requests.
void timeStepRk(const MFloat t, const MFloat dt, const MInt substep=-1)
Time integration using the five-stage fourth-order low-storage Runge-Kutta scheme as described in the...
void saveNodalData(const MString &fileNameBase, const MInt noVars, const std::vector< MString > &varNames, const MFloat *const data) const
Save nodal data to file.
void calcDgTimeDerivative(const MFloat t, const MFloat tStage)
Main routine to calculate the discontinuous Galerkin time derivative. After this method was called,...
void subTimeStepRk(const MFloat dt, const MInt stage, const MInt totalSize, const MFloat *const rhs, MFloat *const variables, MFloat *const timeIntStorage)
Perform one Runge-Kutta substep on the given elements.
void extendNodeVariablesSingleDirection(const MInt extendDir, const MFloat extendOffset, const MFloat extendLimit)
Set nodeVars upstream of given plane to values of elements in plane.
MInt getElementIdAtPoint(const MFloat *point, MBool globalUnique=false)
returns the elementId of a element containing a given point
void initPrefinement()
Set polynomial degree for static p-refinement case.
void initGridMap()
Determine grid-to-grid mapping.
void run()
Main method to run this solver.
MBool prepareRestart(MBool, MBool &) override
Return true if the solver wants to write a restart file.
void finalizeInitSolver() override
Finalization of the solver initialization.
void localToGlobalIds() override
Change local into global ids.
void finishMpiSurfaceExchange()
Finish sending the window-side data and finish receiving the halo-side data for all MPI surfaces.
void initHElements()
Initialize all helements by iterating over all cells and creating an element for each coarse cell in ...
typename maia::CartesianSolver< nDim, DgCartesianSolver > CartesianSolver
void postTimeStep()
Perform post-time-step operations, e.g. advance time, error analysis.
MBool step(const MFloat externalDt=-std::numeric_limits< MFloat >::infinity(), const MBool substep=false)
Advance the solution by one time step.
MBool isAdaptationTimeStep(const MInt timeStep) const
Return if the given timestep is an adaptation timestep.
void resetSolver() override
Reset solver (for load balancing)
void setNumericalProperties()
Reads properties associated with the numerical method.
void calcSurfaceNodeCoordinates(const MInt surfaceId)
void loadRestartFile() override
Load restart file with all information that is necessary to restart the calculations from here.
void balancePre() override
Reinitialize solver prior to setting solution data in DLB.
void initSolverObjects()
Initializes the solver. This must be called before using any of the discretization methods,...
void applyForwardProjection()
void startMpiSurfaceExchange()
Start sending the window-side data and start receiving the halo-side data for all MPI surfaces.
void finalizeMpiExchange()
void updateNodeVariables()
Update all node variables at the surfaces.
MBool isMpiRoot() const
Return true if this is the root rank of the solver MPI communicator.
virtual void initialCondition()
Set the initial condition in all elements.
void saveSolutionFile()
Saves all available data to disk.
void getDefaultWeights(MFloat *weights, std::vector< MString > &names) const
Return the default weights for all load quantities.
void outputInitSummary()
Print initialization summary to user and m_log.
virtual void resetRHS()
Reset the time derivative of the conservative variables to zero.
Class stores precalculated values for interpolation & integration on the reference interval [-1,...
MFloatVector m_LhatFace[2]
void calcNormal(const MFloat *const vertices, MFloat *normal) const
Calculates the normal vector from the geometry element vertices.
void calcCentroid(const MFloat *const vertices, MFloat *centroid) const
Calculate the centroid of the geometry element vertices.
void getVertices(MFloat *vertices) const
Return the vertices of the geometry element.
void open(const MString &filename, const MString &projectName, MInt fileType=0, MPI_Comm mpiComm=MPI_COMM_WORLD, MBool rootOnlyHardwired=false)
Opens a file by passing the parameters to InfoOut_<xyz>FileBuffer::open(...).
This class is a ScratchSpace.
T * getPointer() const
Deprecated: use begin() instead!
pointer p
Deprecated: use [] instead!
MPI_Comm mpiComm() const
Return the MPI communicator used by this solver.
const MInt m_solverId
a unique solver identifier
Class that represents DG element collector.
MInt & polyDeg(const MInt id)
Accessor for polynomial degree.
MInt noNodesXD(const MInt id) const
Accessor for number of nodes XD (const version).
MFloat & invJacobian(const MInt id)
Accessor for inverse jacobian.
MInt & cellId(const MInt id)
Accessor for cell id.
MInt & surfaceIds(const MInt id, const MInt dir)
Accessor for surface ids.
MFloat & rightHandSide(const MInt id)
Accessor for right hand side.
MInt noNodes1D(const MInt id) const
Accessor for number of nodes 1D (const version).
Class that represents DG element collector.
MInt & hrefSurfaceIds(const MInt id, const MInt dir, const MInt pos)
Accessor for h-refined surface ids.
MInt & elementId(const MInt id)
Accessor for element id.
Class that represents DG element collector.
MInt & noNodes1D(const MInt srfcId)
Accessor for number of nodes 1D.
MInt & nghbrElementIds(const MInt srfcId, const MInt side)
Accessor for neighbor element ids.
MFloat & flux(const MInt srfcId)
Accessor for flux.
MInt & fineCellId(const MInt srfcId)
Accessor for fine cell id.
MInt & polyDeg(const MInt srfcId)
Accessor for polynomial degree.
MInt & orientation(const MInt srfcId)
Accessor for orientation.
size_type dim0() const
Return the size of dimension 0.
void clear()
Deletes the data (if internal storage was used) and resets dimensions to zero.
size_type dim1() const
Return the size of dimension 1.
void set(const T &value)
Initializes tensor to constant value.
size_type size() const
Returns the size of the array as product of all five dimensions (i.e. not the actual array size but t...
size_type resize(size_type n0, size_type n1=1, size_type n2=1, size_type n3=1, size_type n4=1)
Deletes the old data structure and creates a new one with the requested dimensions.
size_type dim2() const
Return the size of dimension 2.
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
@ DG_INTEGRATE_GAUSS_LOBATTO
@ DG_TIMEINTEGRATION_CARPENTER_4_5
@ DG_TIMEINTEGRATION_TOULORGEC_3_7
@ DG_TIMEINTEGRATION_TOULORGEC_4_8
@ DG_TIMEINTEGRATION_TOULORGEF_4_8
@ DG_TIMEINTEGRATION_NIEGEMANN_4_13
@ DG_TIMEINTEGRATION_NIEGEMANN_4_14
void mTerm(const MInt errorCode, const MString &location, const MString &message)
const MString const MString & message
MBool approx(const T &, const U &, const T)
MInt ipow(MInt base, MInt exp)
Integer exponent function for non-negative exponents.
void printAllocatedMemory(const MLong oldAllocatedBytes, const MString &solverName, const MPI_Comm comm)
Prints currently allocated memory.
MInt globalDomainId()
Return global domain id.
constexpr MLong IPOW2(MInt x)
constexpr MFloat FPOW2(MInt x)
std::basic_string< char > MString
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Isend
int MPI_Send_init(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Send_init
int MPI_Allgatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, const int recvcounts[], const int displs[], MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgatherv
int MPI_Waitsome(int incount, MPI_Request array_of_requests[], int *outcount, int array_of_indices[], MPI_Status array_of_statuses[], const MString &name)
same as MPI_Waitsome
int MPI_Type_commit(MPI_Datatype *datatype, const MString &name)
same as MPI_Type_commit
int MPI_Irecv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Irecv
int MPI_Recv_init(void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Recv_init
int MPI_Type_free(MPI_Datatype *datatype, const MString &name)
same as MPI_Type_free
int MPI_Cancel(MPI_Request *request, const MString &name)
same as MPI_cancel
int MPI_Wait(MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Wait
int MPI_Waitall(int count, MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Waitall
int MPI_Type_create_hindexed(int count, const int array_of_solverlengths[], const MPI_Aint array_of_displacements[], MPI_Datatype oldtype, MPI_Datatype *newtype, const MString &name)
same as MPI_Type_create_hindexed
int MPI_Get_address(const void *location, MPI_Aint *address, const MString &name)
same as MPI_Get_address
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
int MPI_Request_free(MPI_Request *request, const MString &name)
same as MPI_Request_free
int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm, const MString &name, const MString &varname)
same as MPI_Bcast
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgather
int MPI_Startall(int count, MPI_Request array_of_requests[], const MString &name)
same as MPI_Startall
void calcTrilinearInterpolation(MFloat *point, const MFloat *nodes, const MInt noNodes, const MInt noVars, const MFloat *u, MFloat *const out)
Calculates trilinear interpolation at given point (3D)
void calcLagrangeInterpolatingPolynomials(const MFloat x, const MInt polyDeg, const MFloat *nodes, const MFloat *wBary, MFloat *polynomials)
Calculates the values of the Lagrangian polynomials l_j for a given point x in [-1,...
void calcPolynomialInterpolationMatrix(MInt noNodesIn, const T nodesIn, MInt noNodesOut, const U nodesOut, const V wBary, W vandermonde)
Calculate the polynomial interpolation matrix (Vandermonde) to interpolate from one set of nodes to a...
void calcLinearInterpolationMatrix(const MInt noNodesIn, const T nodesIn, const MInt noNodesOut, const U nodesOut, V vandermonde)
Calculates the linear interpolation matrix (Vandermonde) to interpolate from one set of nodes to anot...
void calcBilinearInterpolation(MFloat *point, const MFloat *nodes, const MInt noNodes, const MInt noVars, const MFloat *u, MFloat *const out)
Calculates bilinear interpolation at given point (2D)
void calcMortarProjectionMatrixP(const MInt sourcePolyDeg, const MFloat *const sourceNodes, const MFloat *const sourceBaryWeights, const MInt targetPolyDeg, const MFloat *const targetNodes, MFloat *const matrix)
T linear(const T a, const T b, const T x)
Linear slope filter.
void calcMortarProjectionMatrixPSBP(const MString sourceOp, const MString targetOp, const MInt sourceNoNodes, const MInt targetNoNodes, MFloat *const f, MFloat *const b)
Reads projection coefficients and stores them.
Namespace for auxiliary functions/classes.
PARALLELIO_DEFAULT_BACKEND ParallelIo
static void init(FactoryType &factory)
Data & addData(const std::string &title_, T data_)
Group & addGroup(const std::string &title_)
Data & addData(const std::string &title_, T data_)
maia::tensor::Tensor< MFloat > MFloatTensor