11#include <unordered_map>
28using namespace std::placeholders;
38 :
maia::CartesianSolver<nDim,
LbSolver<nDim>>(
id, gridProxy_, comm, true)
49 m_geometry(&geometry_),
52 m_noDistributions(
dist) {
56 RECORD_TIMER_START(m_t.solver);
60 auto solverMethodEnum =
string2enum(solverMethod());
65 m_cells.setThermal(
true);
68 m_cells.setThermal(
false);
73 m_cells.setTransport(
true);
76 m_cells.setTransport(
false);
86 m_cells.setNoVariables();
96 m_isEELiquid = Context::getSolverProperty<MBool>(
"EELiquid", m_solverId, AT_, &m_isEELiquid);
98 MInt l_savePrevVars = 0;
99 l_savePrevVars = Context::getBasicProperty<MInt>(
"alphaConvergenceCheck", AT_, &l_savePrevVars);
100 m_cells.setSaveUOtherPhase(
true);
101 m_cells.setSaveVolumeFraction(
true);
102 m_cells.setSavePrevVars(l_savePrevVars > 0);
103 m_cells.setSaveNuT(
true);
104 m_cells.setSaveOldNu(
true);
113 m_EELiquid.restartWithoutAlpha =
false;
114 m_EELiquid.restartWithoutAlpha =
115 Context::getSolverProperty<MBool>(
"LBRestartWithoutAlpha", m_solverId, AT_, &m_EELiquid.restartWithoutAlpha);
119 m_cells.setSaveOldNu(
true);
123 MBool updateAfterPropagation =
false;
125 updateAfterPropagation =
true;
135 updateAfterPropagation =
136 Context::getSolverProperty<MBool>(
"updateAfterPropagation", m_solverId, AT_, &updateAfterPropagation);
137 if(updateAfterPropagation) {
139 TERMM(1,
"For now only implemented for the cumulant collision step!");
146 m_cells.setNoDistributions(m_noDistributions);
151 m_cells.reset(grid().raw().treeb().capacity());
152 m_cells.append(grid().noCells());
155 updateCellCollectorFromGrid();
157 grid().findEqualLevelNeighborsParDiagonal(
false);
160 m_adaptationSinceLastRestart =
false;
161 m_adaptationSinceLastSolution =
false;
171 m_adaptation =
false;
172 m_adaptation = Context::getSolverProperty<MBool>(
"adaptation", m_solverId, AT_, &m_adaptation);
175 mAlloc(this->m_recalcIds, maxNoGridCells(),
"m_recalcIds", -1, AT_);
176 for(
MInt i = 0; i < maxNoGridCells(); i++) {
177 this->m_recalcIds[i] = i;
180 this->m_recalcIds =
nullptr;
192 this->m_singleAdaptation =
true;
193 this->m_singleAdaptation =
194 Context::getSolverProperty<MBool>(
"singleAdaptation", m_solverId, AT_, &this->m_singleAdaptation);
203 m_adaptationInitMethod =
"INIT_FILIPPOVA";
204 m_adaptationInitMethod =
205 Context::getSolverProperty<MString>(
"adaptationInitMethod", m_solverId, AT_, &m_adaptationInitMethod);
208 cerr <<
"adaptation initialization method is: " << m_adaptationInitMethod << endl;
213 m_reinitFileName = this->grid().gridInputFileName().c_str();
214 m_reinitFilePath = outputDir() + m_reinitFileName;
224 noSpecies = Context::getSolverProperty<MInt>(
"noSpecies", m_solverId, AT_, &noSpecies);
235 if(m_isTransport && !m_isThermal) {
236 m_noVariables = nDim + 1 + 2 * m_isTransport;
238 m_noVariables = nDim + 1 + m_isThermal + m_isTransport;
241 m_isRefined = (grid().maxUniformRefinementLevel() < maxLevel());
245 RECORD_TIMER_START(m_t.initSolver);
247 m_log <<
"#########################################################################################################"
250 m_log <<
"## Initializing LB solver "
253 m_log <<
"#########################################################################################################"
268 m_nonBlockingComm =
false;
270 m_nonBlockingComm = Context::getSolverProperty<MBool>(
"nonBlockingComm", m_solverId, AT_, &m_nonBlockingComm);
283 stringstream errorMessage;
284 errorMessage <<
" LbSolver::() Invalid no of distributions : " <<
dist <<
" Exiting...";
285 TERMM(1, errorMessage.str());
288 if(grid().isActive()) {
291 mAlloc(m_activeCellList, grid().noCells(),
"m_activeCellList", AT_);
305 m_solutionOffset = 0;
306 m_solutionOffset = Context::getSolverProperty<MInt>(
"solutionOffset", m_solverId, AT_, &m_solutionOffset);
339 m_initMethod = Context::getSolverProperty<MString>(
"initMethod", m_solverId, AT_);
355 MString interpolationType =
"LINEAR_INTERPOLATION";
357 interpolationType = Context::getSolverProperty<MString>(
"interpolationType", m_solverId, AT_, &interpolationType);
376 m_CouettePoiseuilleRatio = F0;
377 m_CouettePoiseuilleRatio =
378 Context::getSolverProperty<MFloat>(
"CouettePoiseuilleRatio", m_solverId, AT_, &m_CouettePoiseuilleRatio);
392 m_calcTotalPressureGradient = 0;
393 m_calcTotalPressureGradient =
394 Context::getSolverProperty<MInt>(
"calcTotalPressureGradient", m_solverId, AT_, &m_calcTotalPressureGradient);
410 m_densityFluctuations =
false;
411 m_densityFluctuations =
412 Context::getSolverProperty<MBool>(
"densityFluctuations", m_solverId, AT_, &m_densityFluctuations);
428 m_calculateDissipation =
false;
429 m_calculateDissipation =
430 Context::getSolverProperty<MBool>(
"calculateDissipation", m_solverId, AT_, &m_calculateDissipation);
445 m_FftInit = Context::getSolverProperty<MBool>(
"FFTInit", m_solverId, AT_, &m_FftInit);
448 Context::propertyExists(
"fftInterval", 0) ? Context::getSolverProperty<MInt>(
"fftInterval", m_solverId, AT_) : 0;
462 for(
MInt dir = 0; dir < nDim; dir++) {
463 m_arraySize[dir] = 2;
464 m_arraySize[dir] = Context::getSolverProperty<MInt>(
"arraySize", m_solverId, AT_, &m_arraySize[dir], dir);
476 m_noPeakModes = Context::getSolverProperty<MInt>(
"m_noPeakModes", m_solverId, AT_, &m_noPeakModes);
478 if(m_FftInit && m_noPeakModes == 0) {
479 stringstream errorMessage;
480 errorMessage <<
" noPeakModes was not defined for FFTInit, exiting";
481 TERMM(1, errorMessage.str());
496 m_Cs = Context::getSolverProperty<MFloat>(
"smagorinskyConstant", m_solverId, AT_, &m_Cs);
506 m_deltaX = Context::getSolverProperty<MFloat>(
"filterWidth", m_solverId, AT_, &m_deltaX);
514 mAlloc(m_momentumFlux, a_noCells(), nDim * nDim,
"m_momentumFlux", F0, AT_);
516 mAlloc(m_MijMij, a_noCells(),
"m_MijMij", F0, AT_);
517 mAlloc(m_MijLij, a_noCells(),
"m_MijLij", F0, AT_);
537 m_Pr = Context::getSolverProperty<MFloat>(
"Pr", m_solverId, AT_, &m_Pr);
547 m_initTemperatureKelvin = 1.0546;
548 m_initTemperatureKelvin =
549 Context::getSolverProperty<MFloat>(
"initTemperatureKelvin", m_solverId, AT_, &m_initTemperatureKelvin);
561 m_blasiusPos = Context::getSolverProperty<MFloat>(
"blasiusPos", m_solverId, AT_, &m_blasiusPos);
576 m_Pe = Context::getSolverProperty<MFloat>(
"Pe", m_solverId, AT_, &m_Pe);
587 m_initCon = Context::getSolverProperty<MFloat>(
"initCon", m_solverId, AT_, &m_initCon);
600 m_alpha = Context::getSolverProperty<MFloat>(
"alpha", m_solverId, AT_, &m_alpha);
610 m_saveDerivatives =
false;
612 m_saveDerivatives = Context::getSolverProperty<MBool>(
"saveDerivatives", m_solverId, AT_);
624 m_tanhInit = Context::getSolverProperty<MBool>(
"tanhInit", m_solverId, AT_, &m_tanhInit);
635 m_initRe = Context::getSolverProperty<MFloat>(
"initRe", m_solverId, AT_);
643 m_initTime = Context::getSolverProperty<MInt>(
"initTime", m_solverId, AT_);
651 m_initStartTime = Context::getSolverProperty<MInt>(
"initStartTime", m_solverId, AT_);
665 m_Ma = Context::getSolverProperty<MFloat>(
"Ma", m_solverId, AT_, &m_Ma);
677 m_Re = Context::getSolverProperty<MFloat>(
"Re", m_solverId, AT_, &m_Re);
680 m_ReTau = Context::getSolverProperty<MFloat>(
"ReTau", m_solverId, AT_);
682 m_Re = Context::getSolverProperty<MFloat>(
"Re", m_solverId, AT_);
693 m_rho1 = Context::getSolverProperty<MFloat>(
"rho1", m_solverId, AT_, &m_rho1);
703 m_rho2 = Context::getSolverProperty<MFloat>(
"rho2", m_solverId, AT_, &m_rho2);
708 if(grid().isActive()) {
709 mAlloc(m_cellLength, maxLevel() + 1,
"m_cellLength", -F1, AT_);
710 for(
MInt level = 0; level < maxLevel() + 1; level++) {
711 m_cellLength[level] = c_cellLengthAtLevel(level);
733 const MFloat referenceLengthSTL =
734 Context::getSolverProperty<MFloat>(
"referenceLength", m_solverId, AT_, &referenceLengthSTL);
735 const MFloat lengthMaxRefinementLvl = c_cellLengthAtLevel(grid().maxRefinementLevel());
736 m_referenceLengthSTL = referenceLengthSTL;
737 m_referenceLength = referenceLengthSTL / lengthMaxRefinementLvl;
738 ss1 <<
"referenceLength ";
749 m_referenceLength = Context::getSolverProperty<MFloat>(
"referenceLengthLB", m_solverId, AT_);
750 ss1 <<
"referenceLengthLB ";
761 m_referenceLengthSegId = Context::getSolverProperty<MInt>(
"referenceLengthSegId", m_solverId, AT_);
762 m_referenceLength = calculateReferenceLength(m_referenceLengthSegId);
763 ss1 <<
"referenceLengthSegId ";
768 ss2 <<
"One of the following properties must be given: ";
769 ss2 <<
"referenceLength referenceLengthLB referenceLengthSegId";
771 }
else if(chk != 1) {
773 ss2 <<
"Redundant conflicting input properties: ";
778 if(grid().isActive()) {
780 m_domainLength =
FPOW2(maxLevel()) / reductionFactor();
789 const MFloat domainLengthSTL =
790 Context::getSolverProperty<MFloat>(
"domainLength", m_solverId, AT_, &domainLengthSTL);
791 const MFloat lengthMaxRefinementLvl = c_cellLengthAtLevel(grid().maxRefinementLevel());
792 m_domainLength = domainLengthSTL / lengthMaxRefinementLvl;
803 m_domainLength = Context::getSolverProperty<MFloat>(
"domainLengthLB", m_solverId, AT_, &m_domainLength);
807 TERMM(1,
"Redundant conflicting input properties: domainLength, domainLengthLB");
813 mAlloc(m_Fext, m_noDistributions,
"m_Fext", F0, AT_);
815 mAlloc(m_EELiquid.Fg, m_noDistributions,
"m_EELiquid.Fg", F0, AT_);
817 for(
MInt mi = 0; mi < m_noDistributions; mi++) {
818 m_EELiquid.Fg[mi] = F0;
824 TERMM(1,
"Missing property EELiquidGravity or EELiquidGravityAccel!");
833 m_EELiquid.initialAlpha = 0.0;
834 m_EELiquid.initialAlpha =
835 Context::getSolverProperty<MFloat>(
"initialAlpha", m_solverId, AT_, &m_EELiquid.initialAlpha);
844 m_EELiquid.alphaInf = m_EELiquid.initialAlpha;
845 m_EELiquid.alphaInf = Context::getSolverProperty<MFloat>(
"alphaInf", m_solverId, AT_, &m_EELiquid.alphaInf);
847 if(domainId() == 0) cerr <<
"LB Solver EELiquid!" << endl;
857 m_EELiquid.gravity =
false;
858 m_EELiquid.gravity = Context::getSolverProperty<MBool>(
"EELiquidGravity", m_solverId, AT_, &m_EELiquid.gravity);
867 m_initDensityGradient =
false;
868 m_initDensityGradient =
869 Context::getSolverProperty<MBool>(
"initDensityGradient", m_solverId, AT_, &m_initDensityGradient);
871 if(m_EELiquid.gravity) {
873 TERMM(1,
"Missing property EELiquidGravityAccel!");
875 for(
MInt i = 0; i < nDim; i++) {
884 m_EELiquid.gravityAccelM[i] = 0.0;
885 m_EELiquid.gravityAccelM[i] = Context::getSolverProperty<MFloat>(
"EELiquidGravityAccel", m_solverId, AT_, i);
889 m_externalForcing =
false;
891 m_externalForcing =
true;
894 for(
MInt i = 0; i < nDim; i++) {
904 m_volumeAccel[i] += Context::getSolverProperty<MFloat>(
"volumeAcceleration", m_solverId, AT_, i);
916 m_Ga = Context::getSolverProperty<MFloat>(
"Ga", m_solverId, AT_, &m_Ga);
917 const MFloat nu = m_Ma / F1BCS / m_Re * m_referenceLength;
919 m_volumeAccel[1] -= gravity;
933 m_externalForcing = Context::getSolverProperty<MBool>(
"externalForcing", m_solverId, AT_, &m_externalForcing);
935 m_particleMomentumCoupling =
false;
936 m_particleMomentumCoupling =
937 Context::getSolverProperty<MBool>(
"particleMomentumCoupling", m_solverId, AT_, &m_particleMomentumCoupling);
939 m_saveExternalForces =
false;
940 m_saveExternalForces =
941 Context::getSolverProperty<MBool>(
"saveExternalForces", m_solverId, AT_, &m_saveExternalForces);
962 m_velocityControl.dir = -1;
963 m_velocityControl.dir = Context::getSolverProperty<MInt>(
"velocityControl", m_solverId, AT_, &m_velocityControl.dir);
964 if(m_velocityControl.dir < -1 || m_velocityControl.dir > 2) {
965 TERMM(1,
"Invalid Cartesian direction for velocityControl!");
967 if(m_velocityControl.dir >= 0) {
968 m_controlVelocity =
true;
978 m_velocityControl.restart =
false;
979 m_velocityControl.restart =
980 Context::getSolverProperty<MBool>(
"velocityControlRestart", m_solverId, AT_, &m_velocityControl.restart);
990 m_velocityControl.interval = 100;
991 m_velocityControl.interval =
992 Context::getSolverProperty<MInt>(
"velocityControlInterval", m_solverId, AT_, &m_velocityControl.interval);
993 if(m_velocityControl.interval < 1) {
994 TERMM(1,
"Invalid velocityControlInterval!");
1006 m_velocityControl.KT = 1.0;
1007 m_velocityControl.KT =
1008 Context::getSolverProperty<MFloat>(
"velocityControlKT", m_solverId, AT_, &m_velocityControl.KT);
1019 m_velocityControl.KI = 10000.0;
1020 m_velocityControl.KI =
1021 Context::getSolverProperty<MFloat>(
"velocityControlKI", m_solverId, AT_, &m_velocityControl.KI);
1032 m_velocityControl.KD = 10.0;
1033 m_velocityControl.KD =
1034 Context::getSolverProperty<MFloat>(
"velocityControlKD", m_solverId, AT_, &m_velocityControl.KD);
1036 if(m_controlVelocity) {
1037 for(
MInt i = 0; i < nDim; i++) {
1038 m_volumeAccelBase[i] = m_volumeAccel[i];
1049 m_solidLayerExtension =
false;
1050 m_solidLayerExtension = Context::getSolverProperty<MBool>(
"solidLayer", m_solverId, AT_, &m_solidLayerExtension);
1060 m_writeLsData =
false;
1061 m_writeLsData = Context::getSolverProperty<MBool>(
"writeLsData", m_solverId, AT_, &m_writeLsData);
1070 m_useOnlyCollectedLS =
false;
1071 m_useOnlyCollectedLS =
1072 Context::getSolverProperty<MBool>(
"useOnlyCollectedLs", m_solverId, AT_, &m_useOnlyCollectedLS);
1081 m_allowBndryAsG0 =
false;
1082 m_allowBndryAsG0 = Context::getSolverProperty<MBool>(
"allowBndryAsG0", m_solverId, AT_, &m_allowBndryAsG0);
1085 mAlloc(m_rescoordinates, nDim + 1, nDim,
"m_rescoordinates", F0, AT_);
1086 mAlloc(m_residual, nDim + 1,
"m_residual", F0, AT_);
1087 mAlloc(m_tmpResidual, nDim + 1,
"m_tmpResidual", F0, AT_);
1088 mAlloc(m_tmpResidualLvl, nDim + 1,
"m_tmpResidualLvl", 0, AT_);
1089 mAlloc(m_maxResId, nDim + 1,
"m_ResId", 0, AT_);
1095 m_tanhScaleFactor = 1.0 / (1.0 - (tanh(-2.5) + 1.0));
1098 m_nu = m_Ma * LBCS / m_Re * m_referenceLength;
1099 m_omega = 2.0 / (1.0 + 6.0 * m_nu);
1103 m_nonNewtonian =
true;
1104 const MString modelString = Context::getSolverProperty<MString>(
"nonNewtonianModel", m_solverId, AT_);
1107 m_n = Context::getSolverProperty<MFloat>(
"nonNewtonian_n_bndry", m_solverId, AT_);
1109 m_n = Context::getSolverProperty<MFloat>(
"nonNewtonian_n", m_solverId, AT_);
1111 std::stringstream ss;
1112 ss <<
"nonNewtonianModel : " << modelString <<
" is not implemented!" << std::endl;
1117 if(m_isThermal != 0) {
1118 m_kappa = m_nu / m_Pr;
1119 m_omegaT = 2.0 / (1.0 + 6.0 * m_kappa);
1122 if(m_isTransport != 0) {
1123 m_diffusivity = m_nu * (m_Re / m_Pe);
1124 m_omegaD = 2.0 / (1.0 + 6.0 * m_diffusivity);
1128 m_pulsatileFrequency = m_nu * ((m_alpha * m_alpha) / (m_referenceLength * m_referenceLength));
1130 RECORD_TIMER_STOP(m_t.initSolver);
1133 if(m_isRefined || m_adaptation) {
1146 m_correctInterfaceBcCells =
1147 Context::getSolverProperty<MBool>(
"correctInterfaceBcCells", m_solverId, AT_, &m_correctInterfaceBcCells);
1148 if(grid().isActive()) {
1149 treatInterfaceCells();
1154 if(grid().isActive()) {
1155 setActiveCellList();
1159 m_initRestart =
false;
1160 m_initFromCoarse =
false;
1174 m_initRestart = Context::getSolverProperty<MBool>(
"initRestart", m_solverId, AT_);
1188 m_initFromCoarse = Context::getSolverProperty<MBool>(
"initFromCoarse", m_solverId, AT_, &m_initFromCoarse);
1206 m_isInitRun =
false;
1207 m_isInitRun = Context::getSolverProperty<MBool>(
"isInitRun", m_solverId, AT_, &m_isInitRun);
1209 MString lbInterfaceMethod =
"FILIPPOVA";
1210 lbInterfaceMethod = Context::getSolverProperty<MString>(
"interfaceMethod", m_solverId, AT_, &lbInterfaceMethod);
1213 m_log <<
"#########################################################################################################"
1216 m_log <<
"## Methods and flow field init "
1219 m_log <<
"#########################################################################################################"
1223 m_log <<
" + Initializing methods..." << endl;
1224 m_log <<
" - solver method: " << solverMethod() << endl;
1225 m_log <<
" - interface method: " << lbInterfaceMethod << endl << endl;
1229 if(m_isThermal && !m_isTransport) {
1230 m_propagationStepMethod = &LbSolver::propagation_step_thermal;
1231 }
else if(!m_isThermal && m_isTransport) {
1232 m_propagationStepMethod = &LbSolver::propagation_step_transport;
1233 }
else if(m_isThermal && m_isTransport) {
1234 m_propagationStepMethod = &LbSolver::propagation_step_thermaltransport;
1236 m_propagationStepMethod = &LbSolver::propagation_step;
1240 if(m_isThermal && !m_isTransport) {
1241 m_propagationStepMethod = &LbSolver::propagation_step_thermal_vol;
1242 }
else if(!m_isThermal && m_isTransport) {
1243 m_propagationStepMethod = &LbSolver::propagation_step_transport_vol;
1244 }
else if(m_isThermal && m_isTransport) {
1245 m_propagationStepMethod = &LbSolver::propagation_step_thermaltransport_vol;
1247 m_propagationStepMethod = &LbSolver::propagation_step_vol;
1251 TERMM(1,
"Unknown lb interface method: Exiting!");
1255 MInt nonBlockingComm =
static_cast<MInt>(m_nonBlockingComm);
1257 if(nonBlockingComm != 0) {
1258 TERMM(1,
"Due to the solver reconstruction, this exchange method is not available at the moment");
1259 m_exchangeMethod = &LbSolver::exchangeLbNB;
1261 m_exchangeMethod = &LbSolver::exchangeLb;
1264 switch(solverMethodEnum) {
1266 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1267 m_solutionStepMethod = &LbSolver::bgki_collision_step;
1270 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1271 m_solutionStepMethod = &LbSolver::bgki_collision_step_Guo_forcing;
1274 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1275 m_solutionStepMethod = &LbSolver::bgki_init_collision_step;
1278 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1279 m_solutionStepMethod = &LbSolver::bgkc_collision_step;
1280 m_updateMacroscopicLocation =
1282 m_isCompressible =
true;
1285 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1286 m_solutionStepMethod = &LbSolver::bgki_smagorinsky_collision_step;
1289 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1290 m_solutionStepMethod = &LbSolver::bgki_smagorinsky_collision_step2;
1293 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1294 m_solutionStepMethod = &LbSolver::bgki_smago_wall_collision_step;
1297 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1298 m_solutionStepMethod = &LbSolver::bgki_dynamic_smago_collision_step;
1301 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1302 m_solutionStepMethod = &LbSolver::rbgk_dynamic_smago_collision_step;
1305 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1306 m_solutionStepMethod = &LbSolver::bgki_euler_collision_step;
1309 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1310 m_solutionStepMethod = &LbSolver::bgki_thermal_collision_step;
1311 m_isCompressible =
true;
1314 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1315 m_solutionStepMethod = &LbSolver::bgki_innerEnergy_collision_step;
1316 m_isCompressible =
true;
1319 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1320 m_solutionStepMethod = &LbSolver::bgki_totalEnergy_collision_step;
1321 m_isCompressible =
true;
1324 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1325 m_solutionStepMethod = &LbSolver::bgkc_transport_collision_step;
1326 m_isCompressible =
true;
1329 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1330 m_solutionStepMethod = &LbSolver::bgkc_thermal_transport_collision_step;
1331 m_isCompressible =
true;
1334 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1335 m_solutionStepMethod = &LbSolver::bgkc_innerenergy_transport_collision_step;
1336 m_isCompressible =
true;
1339 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1340 m_solutionStepMethod = &LbSolver::bgkc_totalenergy_transport_collision_step;
1341 m_isCompressible =
true;
1344 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1345 m_solutionStepMethod = &LbSolver::rbgk_collision_step;
1348 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1349 m_solutionStepMethod = &LbSolver::rbgk_smagorinsky_collision_step;
1352 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1353 m_solutionStepMethod = &LbSolver::mrt_collision_step;
1354 m_updateMacroscopicLocation =
1356 m_isCompressible =
true;
1359 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1360 m_solutionStepMethod = &LbSolver::mrt2_collision_step;
1363 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1364 m_solutionStepMethod = &LbSolver::clb_collision_step;
1367 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1368 m_solutionStepMethod = &LbSolver::clb_smagorinsky_collision_step;
1371 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1372 m_solutionStepMethod = &LbSolver::mrt_smagorinsky_collision_step;
1375 m_initializeMethod = &LbSolver::initializeLatticeBgk;
1376 m_solutionStepMethod = &LbSolver::cumulant_collision_step;
1377 m_updateMacroscopicLocation =
1379 m_isCompressible =
true;
1382 TERMM(1,
"Unknown LB method: Exiting!");
1391 if(m_restartFile || m_initFromRestartFile) {
1392 m_initMethod =
"FROM_RESTART_FILE";
1393 m_initializeMethod = &LbSolver::loadRestartFile;
1394 if(m_initRestart || m_initFromCoarse) {
1395 m_initializeMethod = &LbSolver::restartInitLb;
1399 getReLambdaAndUrmsInit();
1401 if(m_innerBandWidth !=
nullptr)
mDeallocate(m_innerBandWidth);
1402 mAlloc(m_innerBandWidth, grid().maxRefinementLevel(),
"m_innerBandWidth", F0, AT_);
1403 if(m_outerBandWidth !=
nullptr)
mDeallocate(m_outerBandWidth);
1404 mAlloc(m_outerBandWidth, grid().maxRefinementLevel(),
"m_outerBandWidth", F0, AT_);
1405 if(m_bandWidth !=
nullptr)
mDeallocate(m_bandWidth);
1406 mAlloc(m_bandWidth, grid().maxRefinementLevel(),
"m_bandWidth", 0, AT_);
1420 MFloat distFac[2] = {18.0, 9.0};
1421 for(
MInt i = 0; i < 2; i++) {
1422 distFac[i] = Context::getSolverProperty<MFloat>(
"mbBandWidth", m_solverId, AT_, &distFac[i], i);
1424 m_outerBandWidth[grid().maxRefinementLevel() - 1] = distFac[0] * c_cellLengthAtLevel(grid().maxRefinementLevel());
1425 m_bandWidth[grid().maxRefinementLevel() - 1] = distFac[0];
1426 for(
MInt i = grid().maxRefinementLevel() - 2; i >= 0; i--) {
1427 m_outerBandWidth[i] = m_outerBandWidth[i + 1] + (distFac[1] * c_cellLengthAtLevel(i + 1));
1428 m_bandWidth[i] = (m_bandWidth[i + 1] / 2) + 1 + distFac[1];
1430 for(
MInt i = 0; i < grid().maxRefinementLevel(); i++) {
1431 m_innerBandWidth[i] = -m_outerBandWidth[i];
1432 m_log <<
"bandwidth level " << i <<
": " << m_innerBandWidth[i] <<
" " << m_outerBandWidth[i] << endl;
1437 if(domainId() == 0) {
1439 m_resFileName =
"maia_res";
1440 m_resFileName +=
".txt";
1441 mRes.
open(m_resFileName, std::ofstream::app);
1444 mRes <<
"This is the residual file of the LB solver." << std::endl;
1446 mRes <<
"#-----------------------------------------------------------------" << std::endl;
1447 mRes <<
"#---SOLVER INFORMATION " << std::endl;
1449 mRes <<
"Number of ranks: " << noDomains() << std::endl;
1450 mRes <<
"Used solver method: " << solverMethod() << std::endl;
1451 mRes <<
"Used initialization method: " << m_initMethod << std::endl;
1452 mRes <<
"Number dimensions: " << nDim << std::endl;
1453 mRes <<
"Number of distributions: " << m_noDistributions << std::endl;
1454 mRes <<
"Minimal spatial refinement level: " << minLevel() << std::endl;
1455 mRes <<
"Maximal spatial refinement level: " << maxLevel() << std::endl;
1457 mRes <<
"#-----------------------------------------------------------------" << std::endl;
1458 mRes <<
"#---RESIDUAL INFORMATION " << std::endl;
1460 mRes <<
"Residual intervall: " << m_residualInterval << std::endl;
1462 mRes <<
"#-----------------------------------------------------------------" << std::endl;
1463 mRes <<
"#---FLOW PARAMETERS " << std::endl;
1465 mRes <<
"Mach number in the simulation: " << m_Ma << std::endl;
1466 mRes <<
"Reynolds number in the simulation: " << m_Re << std::endl;
1468 mRes <<
"#-----------------------------------------------------------------" << std::endl;
1469 mRes << std::endl << std::endl;
1470 mRes <<
"Calculated resudial for each interval: " << std::endl;
1498 if(m_interfaceChildren.size() != 0) {
1499 m_interfaceChildren.clear();
1501 if(m_interfaceParents.size() != 0) {
1502 m_interfaceParents.clear();
1506 RECORD_TIMER_STOP(m_t.solver);
1519 NEW_TIMER_GROUP(tg_solver,
"LB Solver (solverId=" + std::to_string(m_solverId) +
")");
1520 NEW_TIMER_NOCREATE(m_t.solver,
"complete solver", tg_solver);
1522 NEW_SUB_TIMER_NOCREATE(m_t.initSolver,
"init solver", m_t.solver);
1523 NEW_SUB_TIMER_NOCREATE(m_t.solutionStep,
"SolutionStep", m_t.solver);
1524 NEW_SUB_TIMER_NOCREATE(m_t.collision,
"Collision", m_t.solutionStep);
1525 NEW_SUB_TIMER_NOCREATE(m_t.collisionBC,
"CollisionBC", m_t.solutionStep);
1526 NEW_SUB_TIMER_NOCREATE(m_t.propagation,
"Propagation", m_t.solutionStep);
1527 NEW_SUB_TIMER_NOCREATE(m_t.propagationBC,
"PropagationBC", m_t.solutionStep);
1528 NEW_SUB_TIMER_NOCREATE(m_t.exchange,
"Exchange", m_t.solutionStep);
1529 NEW_SUB_TIMER_NOCREATE(m_t.packing,
"Packing", m_t.exchange);
1530 NEW_SUB_TIMER_NOCREATE(m_t.unpacking,
"Unpacking", m_t.exchange);
1531 NEW_SUB_TIMER_NOCREATE(m_t.communication,
"Communication", m_t.exchange);
1532 NEW_SUB_TIMER_NOCREATE(m_t.residual,
"Residuum", m_t.solutionStep);
1533 NEW_SUB_TIMER_NOCREATE(m_t.srcTerms,
"SourceTerms", m_t.solutionStep);
1535 NEW_SUB_TIMER_NOCREATE(m_t.findG0Cells,
"Find G0 cells", m_t.solver);
1536 NEW_SUB_TIMER_NOCREATE(m_t.resetListsMb,
"Reset lists MB", m_t.findG0Cells);
1537 NEW_SUB_TIMER_NOCREATE(m_t.findG0Candidates,
"Find G0 candidates", m_t.findG0Cells);
1538 NEW_SUB_TIMER_NOCREATE(m_t.geomNodal,
"Geometry intersection compute Nodal", m_t.findG0Cells);
1539 NEW_SUB_TIMER_NOCREATE(m_t.geomExchange,
"Geometry intersection exchange", m_t.findG0Cells);
1540 NEW_SUB_TIMER_NOCREATE(m_t.calcNodalValues,
"Calc nodal values", m_t.findG0Cells);
1542 NEW_SUB_TIMER_NOCREATE(m_t.prepComm,
"prepare communication", m_t.solver);
1544 NEW_SUB_TIMER_NOCREATE(m_t.fft,
"FFT", m_t.solver);
1554 if(!grid().isActive())
return;
1558 m_timerType = Context::getSolverProperty<MString>(
"timerType", m_solverId, AT_, &m_timerType);
1561 std::vector<MInt> timerIds_;
1562 timerIds_.reserve(17);
1563 timerIds_.emplace_back(m_t.solver);
1564 timerIds_.emplace_back(m_t.initSolver);
1565 timerIds_.emplace_back(m_t.solutionStep);
1566 timerIds_.emplace_back(m_t.collision);
1567 timerIds_.emplace_back(m_t.collisionBC);
1568 timerIds_.emplace_back(m_t.propagation);
1569 timerIds_.emplace_back(m_t.propagationBC);
1570 timerIds_.emplace_back(m_t.srcTerms);
1571 timerIds_.emplace_back(m_t.exchange);
1572 timerIds_.emplace_back(m_t.residual);
1573 timerIds_.emplace_back(m_t.findG0Cells);
1574 timerIds_.emplace_back(m_t.resetListsMb);
1575 timerIds_.emplace_back(m_t.findG0Candidates);
1576 timerIds_.emplace_back(m_t.geomNodal);
1577 timerIds_.emplace_back(m_t.geomExchange);
1578 timerIds_.emplace_back(m_t.calcNodalValues);
1579 timerIds_.emplace_back(m_t.prepComm);
1580 timerIds_.emplace_back(m_t.fft);
1581 const MInt noTimers = timerIds_.size();
1584 std::vector<MFloat> timerValues_;
1585 timerValues_.reserve(noTimers);
1586 for(
MInt i = 0; i < noTimers; i++) {
1587 timerValues_.emplace_back(RETURN_TIMER_TIME(timerIds_[i]));
1591 if(m_timerType ==
"average") {
1593 AT_,
"MPI_IN_PLACE",
"timerValues_");
1596 AT_,
"MPI_IN_PLACE",
"timerValues_");
1600 if(m_timerType ==
"average") {
1601 const MInt noDomains_ = noDomains();
1602 for(
MInt i = 0; i < noTimers; i++) {
1603 const MFloat meanValue = timerValues_[i] / noDomains_;
1604 SET_RECORD(timerIds_[i], meanValue);
1607 for(
MInt i = 0; i < noTimers; i++) {
1608 SET_RECORD(timerIds_[i], timerValues_[i]);
1615 std::map<MString, MInt> scalingVars;
1617 scalingVars[
"noCells"] = a_noCells();
1618 scalingVars[
"noInternalCells"] = noInternalCells();
1619 scalingVars[
"noHaloCells"] = a_noCells() - grid().noInternalCells();
1620 scalingVars[
"noActiveCells"] = m_currentMaxNoCells;
1622 std::vector<MInt> recvBuffer;
1623 recvBuffer.resize(noDomains());
1624 for(
const auto& [name, value] : scalingVars) {
1626 0, mpiComm(), AT_,
"send",
"recv");
1631 avgValue = std::accumulate(recvBuffer.begin(), recvBuffer.end(), avgValue);
1632 avgValue /= noDomains();
1633 std::cout <<
"Average " + name +
" per Domain are: " << avgValue << std::endl;
1636 const MInt maxValue = *std::max_element(recvBuffer.begin(), recvBuffer.end());
1637 std::cout <<
"Max " + name +
" per Domain are: " << maxValue << std::endl;
1640 const MInt minValue = *std::min_element(recvBuffer.begin(), recvBuffer.end());
1641 std::cout <<
"Min " + name +
" per Domain are: " << minValue << std::endl;
1653 for(
MInt i = 0; i < a_noCells(); i++) {
1654 m_cells.level(i) = c_level(i);
1655 ASSERT(a_level(i) = c_level(i),
"");
1656 a_isHalo(i) =
false;
1657 a_isWindow(i) =
false;
1659 for(
MInt i = 0; i < noNeighborDomains(); i++) {
1660 for(
MInt c = 0; c < noHaloCells(i); c++) {
1661 a_isHalo(haloCell(i, c)) =
true;
1664 for(
MInt i = 0; i < noNeighborDomains(); i++) {
1665 for(
MInt c = 0; c < noWindowCells(i); c++) {
1666 a_isWindow(windowCell(i, c)) =
true;
1688 RECORD_TIMER_START(m_t.prepComm);
1691 m_log <<
"#########################################################################################################"
1694 m_log <<
"## Communication "
1697 m_log <<
"#########################################################################################################"
1712 m_reducedComm =
false;
1713 m_reducedComm = Context::getSolverProperty<MBool>(
"reducedComm", m_solverId, AT_, &m_reducedComm);
1717 m_gatherMethod = &LbSolver::gatherReduced;
1718 m_scatterMethod = &LbSolver::scatterReduced;
1719 m_sendMethod = &LbSolver::sendReduced;
1720 m_receiveMethod = &LbSolver::receiveReduced;
1721 prepareCommunicationReduced();
1724 m_gatherMethod = &LbSolver::gatherNormal;
1725 m_scatterMethod = &LbSolver::scatterNormal;
1726 m_sendMethod = &LbSolver::sendNormal;
1727 m_receiveMethod = &LbSolver::receiveNormal;
1728 prepareCommunicationNormal();
1730 m_log <<
" lbsolver: domain " << domainId() <<
", has " << a_noCells() <<
" cells" << endl;
1732 RECORD_TIMER_STOP(m_t.prepComm);
1761 m_log <<
" + Preparing normal communication ..." << endl << endl;
1764 if(m_isRefined || m_isEELiquid) {
1766 if(m_isThermal && !m_isTransport) {
1767 m_noElementsTransfer = 4;
1768 mAlloc(m_dataBlockSizes, m_noElementsTransfer,
"m_noDataBlockSizes", 0, AT_);
1769 mAlloc(m_baseAddresses, m_noElementsTransfer,
"m_baseAddresses", AT_);
1771 m_dataBlockSizes[0] = m_noDistributions;
1772 m_dataBlockSizes[1] = m_noDistributions;
1773 m_dataBlockSizes[2] = nDim + 2;
1774 m_dataBlockSizes[3] = nDim + 2;
1776 m_baseAddresses[0] = &a_distribution(0, 0);
1777 m_baseAddresses[1] = &a_distributionThermal(0, 0);
1778 m_baseAddresses[2] = &a_variable(0, 0);
1779 m_baseAddresses[3] = &a_oldVariable(0, 0);
1780 }
else if(m_isTransport && !m_isThermal) {
1781 m_noElementsTransfer = 4;
1782 mAlloc(m_dataBlockSizes, m_noElementsTransfer,
"m_noDataBlockSizes", 0, AT_);
1783 mAlloc(m_baseAddresses, m_noElementsTransfer,
"m_baseAddresses", AT_);
1785 m_dataBlockSizes[0] = m_noDistributions;
1786 m_dataBlockSizes[1] = m_noDistributions;
1787 m_dataBlockSizes[2] = nDim + 2;
1788 m_dataBlockSizes[3] = nDim + 2;
1790 m_baseAddresses[0] = &a_distribution(0, 0);
1791 m_baseAddresses[1] = &a_distributionTransport(0, 0);
1792 m_baseAddresses[2] = &a_variable(0, 0);
1793 m_baseAddresses[3] = &a_oldVariable(0, 0);
1794 }
else if(m_isThermal && m_isTransport) {
1795 m_noElementsTransfer = 5;
1796 mAlloc(m_dataBlockSizes, m_noElementsTransfer,
"m_noDataBlockSizes", 0, AT_);
1797 mAlloc(m_baseAddresses, m_noElementsTransfer,
"m_baseAddresses", AT_);
1799 m_dataBlockSizes[0] = m_noDistributions;
1800 m_dataBlockSizes[1] = m_noDistributions;
1801 m_dataBlockSizes[2] = m_noDistributions;
1802 m_dataBlockSizes[3] = nDim + 3;
1803 m_dataBlockSizes[4] = nDim + 3;
1805 m_baseAddresses[0] = &a_distribution(0, 0);
1806 m_baseAddresses[1] = &a_distributionThermal(0, 0);
1807 m_baseAddresses[2] = &a_distributionTransport(0, 0);
1808 m_baseAddresses[3] = &a_variable(0, 0);
1809 m_baseAddresses[4] = &a_oldVariable(0, 0);
1813 m_noElementsTransfer = 3;
1814 mAlloc(m_dataBlockSizes, m_noElementsTransfer,
"m_noDataBlockSizes", 0, AT_);
1815 mAlloc(m_baseAddresses, m_noElementsTransfer,
"m_baseAddresses", AT_);
1817 m_dataBlockSizes[0] = m_noDistributions;
1818 m_dataBlockSizes[1] = nDim + 1;
1819 m_dataBlockSizes[2] = nDim + 1;
1821 m_baseAddresses[0] = &a_distribution(0, 0);
1822 m_baseAddresses[1] = &a_variable(0, 0);
1823 m_baseAddresses[2] = &a_oldVariable(0, 0);
1827 if(m_isThermal && !m_isTransport) {
1828 m_noElementsTransfer = 2;
1829 mAlloc(m_dataBlockSizes, m_noElementsTransfer,
"m_noDataBlockSizes", 0, AT_);
1830 mAlloc(m_baseAddresses, m_noElementsTransfer,
"m_baseAddresses", AT_);
1832 m_dataBlockSizes[0] = m_noDistributions;
1833 m_dataBlockSizes[1] = m_noDistributions;
1835 m_baseAddresses[0] = &a_distribution(0, 0);
1836 m_baseAddresses[1] = &a_distributionThermal(0, 0);
1837 }
else if(m_isTransport && !m_isThermal) {
1838 m_noElementsTransfer = 2;
1839 mAlloc(m_dataBlockSizes, m_noElementsTransfer,
"m_noDataBlockSizes", 0, AT_);
1840 mAlloc(m_baseAddresses, m_noElementsTransfer,
"m_baseAddresses", AT_);
1842 m_dataBlockSizes[0] = m_noDistributions;
1843 m_dataBlockSizes[1] = m_noDistributions;
1845 m_baseAddresses[0] = &a_distribution(0, 0);
1846 m_baseAddresses[1] = &a_distributionTransport(0, 0);
1848 }
else if(m_isTransport && m_isThermal) {
1849 m_noElementsTransfer = 3;
1850 mAlloc(m_dataBlockSizes, m_noElementsTransfer,
"m_noDataBlockSizes", 0, AT_);
1851 mAlloc(m_baseAddresses, m_noElementsTransfer,
"m_baseAddresses", AT_);
1853 m_dataBlockSizes[0] = m_noDistributions;
1854 m_dataBlockSizes[1] = m_noDistributions;
1855 m_dataBlockSizes[2] = m_noDistributions;
1857 m_baseAddresses[0] = &a_distribution(0, 0);
1858 m_baseAddresses[1] = &a_distributionThermal(0, 0);
1859 m_baseAddresses[2] = &a_distributionTransport(0, 0);
1865 m_noElementsTransfer = 2;
1866 mAlloc(m_dataBlockSizes, m_noElementsTransfer,
"m_noDataBlockSizes", 0, AT_);
1867 mAlloc(m_baseAddresses, m_noElementsTransfer,
"m_baseAddresses", AT_);
1869 m_dataBlockSizes[0] = m_noDistributions;
1870 m_dataBlockSizes[1] = nDim + 1;
1872 m_baseAddresses[0] = &a_distribution(0, 0);
1873 m_baseAddresses[1] = &a_variable(0, 0);
1878 m_dataBlockSizeTotal = 0;
1879 for(
MInt i = 0; i < m_noElementsTransfer; i++)
1880 m_dataBlockSizeTotal += m_dataBlockSizes[i];
1883 if(noNeighborDomains() > 0) {
1884 mAlloc(m_nghbrOffsetsWindow, noNeighborDomains(), m_noElementsTransfer,
"m_nghbrOffsetsWindow", 0, AT_);
1885 mAlloc(m_nghbrOffsetsHalo, noNeighborDomains(), m_noElementsTransfer,
"m_nghbrOffsetsHalo", 0, AT_);
1887 for(
MInt n = 0; n < noNeighborDomains(); n++) {
1888 for(
MInt v = 1; v < m_noElementsTransfer; v++) {
1889 m_nghbrOffsetsWindow[n][v] = m_nghbrOffsetsWindow[n][v - 1] + noWindowCells(n) * m_dataBlockSizes[v - 1];
1890 m_nghbrOffsetsHalo[n][v] = m_nghbrOffsetsHalo[n][v - 1] + noHaloCells(n) * m_dataBlockSizes[v - 1];
1894 if(noNeighborDomains() > 0) {
1898 for(
MInt d = 0; d < noNeighborDomains(); d++) {
1899 haloCellsCnt[d] = noHaloCells(d);
1900 windowCellsCnt[d] = noWindowCells(d);
1902 mAlloc(m_sendBuffers, noNeighborDomains(), &windowCellsCnt[0], m_dataBlockSizeTotal,
"m_sendBuffers", AT_);
1903 mAlloc(m_receiveBuffers, noNeighborDomains(), &haloCellsCnt[0], m_dataBlockSizeTotal,
"m_receiveBuffers", AT_);
1906 if(!m_nonBlockingComm) {
1907 mAlloc(mpi_request, noNeighborDomains(),
"mpi_request", AT_);
1911 mAlloc(mpi_requestS, noNeighborDomains(),
"mpi_requestS", AT_);
1912 mAlloc(mpi_requestR, noNeighborDomains(),
"mpi_requestR", AT_);
1917 printCommunicationMethod();
1936 m_log <<
" + Preparing reduced communication ..." << endl << endl;
1939 if(noNeighborDomains() > 0) {
1940 if(!m_nonBlockingComm) {
1941 mAlloc(mpi_request, noNeighborDomains(),
"mpi_request", AT_);
1943 mAlloc(mpi_requestS, noNeighborDomains(),
"mpi_requestS", AT_);
1944 mAlloc(mpi_requestR, noNeighborDomains(),
"mpi_requestR", AT_);
1948 ScratchSpace<MInt> noWindowsPerDomain(noNeighborDomains(), AT_,
"noWindowsPerDomain");
1949 for(
MInt n = 0; n < noNeighborDomains(); n++) {
1950 noHalosPerDomain[n] = noHaloCells(n) * (m_noDistributions - 1) + 1;
1951 noWindowsPerDomain[n] = noWindowCells(n) * (m_noDistributions - 1) + 1;
1954 mAlloc(m_noWindowDistDataPerDomain, noNeighborDomains(),
"m_noWindowDistDataPerDomain", 0, AT_);
1955 mAlloc(m_noHaloDistDataPerDomain, noNeighborDomains(),
"m_noHaloDistDataPerDomain", 0, AT_);
1957 mAlloc(m_windowDistsForExchange, noNeighborDomains(), &noWindowsPerDomain[0],
"m_windowDistsForExchange", AT_);
1958 mAlloc(m_haloDistsForExchange, noNeighborDomains(), &noHalosPerDomain[0],
"m_haloDistsForExchange", AT_);
1961 mAlloc(m_needsFurtherExchange, a_noCells(),
"needsFurtherExchange", 0, AT_);
1962 markCellsForAdditionalComm();
1966 m_log <<
" - counting PPDFs per neighbor" << endl;
1968 for(
MInt n = 0; n < noNeighborDomains(); n++) {
1969 MInt haloStart = domainOffset(neighborDomain(n));
1970 MInt haloEnd = domainOffset(neighborDomain(n) + 1);
1973 for(
MInt j = 0; j < noWindowCells(n); j++) {
1974 if(m_needsFurtherExchange[windowCell(n, j)] > 0) {
1975 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
1976 m_noWindowDistDataPerDomain[n]++;
1977 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] = cnt;
1981 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
1982 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] = -1;
1983 MInt ngh = c_neighborId(windowCell(n, j), d);
1984 if(ngh != -1 && a_isHalo(ngh)) {
1985 MInt global = c_globalId(ngh);
1986 if(global >= haloStart && global < haloEnd) {
1987 m_noWindowDistDataPerDomain[n]++;
1988 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] = cnt;
1996 for(
MInt j = 0; j < noHaloCells(n); j++) {
1997 if(m_needsFurtherExchange[haloCell(n, j)] > 0) {
1998 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
1999 m_noHaloDistDataPerDomain[n]++;
2000 m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] = cnt;
2004 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
2005 m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] = -1;
2006 MInt ngh = c_neighborId(haloCell(n, j), d);
2007 if(ngh != -1 && a_isWindow(ngh)) {
2008 m_noHaloDistDataPerDomain[n]++;
2009 m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] = cnt;
2016 m_log <<
" * neighbor " << neighborDomain(n) <<
":" << endl;
2017 m_log <<
" = windows :" << m_noWindowDistDataPerDomain[n] << endl;
2018 m_log <<
" = halos :" << m_noHaloDistDataPerDomain[n] << endl;
2022 if(m_isThermal && !m_isTransport) {
2023 m_noDistsTransfer = 2;
2024 m_noVarsTransfer = 2;
2025 m_noElementsTransfer = m_noDistsTransfer + m_noVarsTransfer;
2026 }
else if(m_isTransport && !m_isThermal) {
2027 m_noDistsTransfer = 2;
2028 m_noVarsTransfer = 2;
2029 m_noElementsTransfer = m_noDistsTransfer + m_noVarsTransfer;
2030 }
else if(m_isTransport && m_isThermal) {
2031 m_noDistsTransfer = 3;
2032 m_noVarsTransfer = 2;
2033 m_noElementsTransfer = m_noDistsTransfer + m_noVarsTransfer;
2035 m_noDistsTransfer = 1;
2036 m_noVarsTransfer = 2;
2037 m_noElementsTransfer = m_noDistsTransfer + m_noVarsTransfer;
2041 m_log <<
" - counting PPDFs per neighbor" << endl;
2043 for(
MInt n = 0; n < noNeighborDomains(); n++) {
2044 MInt haloStart = domainOffset(neighborDomain(n));
2045 MInt haloEnd = domainOffset(neighborDomain(n) + 1);
2048 for(
MInt j = 0; j < noWindowCells(n); j++) {
2049 if(!c_isLeafCell(windowCell(n, j))) {
2050 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
2051 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] = -1;
2054 if(m_needsFurtherExchange[windowCell(n, j)] > 0) {
2055 for(
MInt d = 0; d < m_noDistributions; d++) {
2056 m_noWindowDistDataPerDomain[n]++;
2057 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] = cnt;
2061 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
2062 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] = -1;
2063 MInt ngh = c_neighborId(windowCell(n, j), d);
2064 if(ngh != -1 && a_isHalo(ngh)) {
2065 MInt global = c_globalId(ngh);
2066 if(global >= haloStart && global < haloEnd) {
2067 m_noWindowDistDataPerDomain[n]++;
2068 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] = cnt;
2077 for(
MInt j = 0; j < noHaloCells(n); j++) {
2078 if(!c_isLeafCell(haloCell(n, j))) {
2079 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
2080 m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] = -1;
2083 if(m_needsFurtherExchange[haloCell(n, j)] > 0) {
2084 for(
MInt d = 0; d < m_noDistributions; d++) {
2085 m_noHaloDistDataPerDomain[n]++;
2086 m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] = cnt;
2090 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
2091 m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] = -1;
2092 MInt ngh = c_neighborId(haloCell(n, j), d);
2093 if(ngh != -1 && a_isWindow(ngh)) {
2094 m_noHaloDistDataPerDomain[n]++;
2095 m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] = cnt;
2103 m_log <<
" * neighbor " << neighborDomain(n) <<
":" << endl;
2104 m_log <<
" = windows :" << m_noWindowDistDataPerDomain[n] << endl;
2105 m_log <<
" = halos :" << m_noHaloDistDataPerDomain[n] << endl;
2109 if(m_isThermal && !m_isTransport) {
2110 m_noDistsTransfer = 2;
2111 m_noVarsTransfer = 0;
2112 m_noElementsTransfer = m_noDistsTransfer + m_noVarsTransfer;
2113 }
else if(m_isTransport && !m_isThermal) {
2114 m_noDistsTransfer = 2;
2115 m_noVarsTransfer = 0;
2116 m_noElementsTransfer = m_noDistsTransfer + m_noVarsTransfer;
2117 }
else if(m_isTransport && m_isThermal) {
2118 m_noDistsTransfer = 3;
2119 m_noVarsTransfer = 0;
2120 m_noElementsTransfer = m_noDistsTransfer + m_noVarsTransfer;
2122 m_noDistsTransfer = 1;
2123 m_noVarsTransfer = 1;
2124 m_noElementsTransfer = m_noDistsTransfer + m_noVarsTransfer;
2129 if(noNeighborDomains() > 0) {
2133 for(
MInt n = 0; n < noNeighborDomains(); n++) {
2134 windowDataPerDomain(n) =
2135 m_noWindowDistDataPerDomain[n] * m_noDistsTransfer + noWindowCells(n) * m_noVariables * m_noVarsTransfer;
2136 haloDataPerDomain(n) =
2137 m_noHaloDistDataPerDomain[n] * m_noDistsTransfer + noHaloCells(n) * m_noVariables * m_noVarsTransfer;
2140 mAlloc(m_sendBuffers, noNeighborDomains(), &windowDataPerDomain[0],
"m_sendBuffers", AT_);
2141 mAlloc(m_receiveBuffers, noNeighborDomains(), &haloDataPerDomain[0],
"m_receiveBuffers", AT_);
2145 printCommunicationMethod();
2164 MInt periodicSimulation = 0;
2165 for(
MInt i = 0; i < a_noCells(); i++) {
2166 if(grid().isPeriodic(i)) {
2167 periodicSimulation = 1;
2171 MPI_Allreduce(MPI_IN_PLACE, &periodicSimulation, 1, MPI_INT, MPI_MAX, mpiComm(), AT_,
"MPI_IN_PLACE",
2172 "periodicSimulation");
2174 if(periodicSimulation) {
2175 MInt noPeriodicHalos = 0;
2176 MIntScratchSpace noPeriodicHalosPerDomain(noNeighborDomains(), AT_,
"noDomains");
2177 std::vector<MInt> myGlobIdsPeriHalos;
2178 for(
MInt i = 0; i < a_noCells(); i++) {
2179 if(!a_isHalo(i))
continue;
2180 if(grid().isPeriodic(i)) {
2181 myGlobIdsPeriHalos.push_back(c_globalId(i));
2186 for(
MInt i = 0; i < noNeighborDomains(); i++) {
2187 MPI_Issend(&noPeriodicHalos, 1, MPI_INT, neighborDomain(i), 0, mpiComm(), &mpi_request[i], AT_,
2191 for(
MInt i = 0; i < noNeighborDomains(); i++) {
2192 MPI_Recv(&noPeriodicHalosPerDomain(i), 1, MPI_INT, neighborDomain(i), 0, mpiComm(), &status, AT_,
2193 "noPeriodicHalosPerDomain(i)");
2195 for(
MInt i = 0; i < noNeighborDomains(); i++) {
2196 MPI_Wait(&mpi_request[i], &status, AT_);
2198 MInt** allGlobIdsPeriHalos{};
2199 mAlloc(allGlobIdsPeriHalos, noNeighborDomains(), &noPeriodicHalosPerDomain[0],
"allGlobalIdsPeriHalos", AT_);
2201 for(
MInt i = 0; i < noNeighborDomains(); i++) {
2202 MPI_Issend(myGlobIdsPeriHalos.data(), noPeriodicHalos, MPI_INT, neighborDomain(i), 0, mpiComm(), &mpi_request[i],
2203 AT_,
"noPeriodicHalos");
2205 for(
MInt i = 0; i < noNeighborDomains(); i++) {
2206 MPI_Recv(allGlobIdsPeriHalos[i], noPeriodicHalosPerDomain(i), MPI_INT, neighborDomain(i), 0, mpiComm(), &status,
2207 AT_,
"noPeriodicHalosPerDomain(i)");
2209 for(
MInt i = 0; i < noNeighborDomains(); i++) {
2210 MPI_Wait(&mpi_request[i], &status, AT_);
2213 for(
MInt i = 0; i < a_noCells(); i++) {
2214 if(a_isHalo(i))
continue;
2215 MInt globalId = c_globalId(i);
2216 MBool found =
false;
2217 for(
MInt d = 0; d < noNeighborDomains(); d++) {
2218 for(
MInt halo = 0; halo < noPeriodicHalosPerDomain(d); halo++) {
2219 if(allGlobIdsPeriHalos[d][halo] == globalId) {
2220 m_needsFurtherExchange[i] = 1;
2230 for(
MInt i = 0; i < a_noCells(); i++) {
2231 if(a_isHalo(i))
continue;
2232 if(a_isInterfaceParent(i)) m_needsFurtherExchange[i] = 1;
2233 if(a_isInterfaceChild(i)) m_needsFurtherExchange[i] = 1;
2234 if(a_isBndryCell(i)) m_needsFurtherExchange[i] = 1;
2235 if(c_parentId(i) > 0) {
2236 if(a_isInterfaceParent(c_parentId(i))) m_needsFurtherExchange[i] = 1;
2238 for(
MInt n = 0; n < IPOW3[nDim] - 1; n++) {
2239 if(c_neighborId(i, n) > -1) {
2240 if(a_isActive(i) != a_isActive(c_neighborId(i, n))) {
2241 m_needsFurtherExchange[i] = 1;
2250 for(
MInt i = 0; i < a_noCells(); i++) {
2251 if(a_isHalo(i))
continue;
2252 if(m_needsFurtherExchange[i])
continue;
2253 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
2254 MInt nghbr = c_neighborId(i, d);
2256 if(m_needsFurtherExchange[nghbr] == 1) m_needsFurtherExchange[i] = 2;
2263 for(
MInt i = 0; i < a_noCells(); i++) {
2264 if(a_isHalo(i))
continue;
2265 if(m_needsFurtherExchange[i])
continue;
2266 for(
MInt d = 0; d < m_noDistributions - 1; d++) {
2267 MInt nghbr = c_neighborId(i, d);
2269 if(m_needsFurtherExchange[nghbr] == 2) m_needsFurtherExchange[i] = 3;
2287 if(!m_nonBlockingComm)
2288 m_log <<
" + Communication method: blocking " << (m_reducedComm ?
" (reduced)" :
" (normal)") << endl << endl;
2290 m_log <<
" + Communication method: non-blocking" << (m_reducedComm ?
" (reduced)" :
" (normal)") << endl << endl;
2292 m_log <<
" + Elements to be transferred: " << m_noElementsTransfer <<
" ( ";
2295 m_log <<
"PVs, oldPVs, PPDFs";
2297 m_log <<
", thPPDFs )" << endl;
2299 m_log <<
" )" << endl;
2303 m_log <<
", thPPDFs )" << endl;
2305 m_log <<
" )" << endl;
2308 for(
MInt i = 0; i < m_noElementsTransfer; i++)
2309 m_log <<
" - size [" << i <<
"] : " << m_dataBlockSizes[i] << endl;
2311 if(!m_reducedComm)
m_log <<
" - total size: " << m_dataBlockSizeTotal << endl << endl;
2313 m_log <<
" + Neighboring domains are:" << endl;
2315 for(
MInt i = 0; i < noNeighborDomains(); i++) {
2316 m_log <<
" - domain id " << neighborDomain(i) <<
":" << endl;
2317 m_log <<
" * window: " << noWindowCells(i) << endl;
2318 m_log <<
" * halo: " << noHaloCells(i) << endl;
2319 m_log <<
" * window data offsets: ";
2321 if(!m_reducedComm) {
2322 for(
MInt v = 0; v < m_noElementsTransfer; v++)
2323 m_log << m_nghbrOffsetsWindow[i][v] <<
" ";
2326 m_log <<
" * halo data offsets: ";
2327 for(
MInt v = 0; v < m_noElementsTransfer; v++)
2328 m_log << m_nghbrOffsetsHalo[i][v] <<
" ";
2332 if(!m_reducedComm) {
2333 m_log <<
" * total window buffer: " << (noWindowCells(i) * m_dataBlockSizeTotal *
sizeof(
MFloat)) <<
" bytes"
2335 m_log <<
" * total halo buffer : " << (noHaloCells(i) * m_dataBlockSizeTotal *
sizeof(
MFloat)) <<
" bytes"
2338 m_log <<
" * total window buffer: " << m_noWindowDistDataPerDomain[i] *
sizeof(
MFloat) <<
" bytes" << endl;
2339 m_log <<
" * total halo buffer : " << m_noHaloDistDataPerDomain[i] *
sizeof(
MFloat) <<
" bytes" << endl;
2343 m_log << endl << endl;
2361 for(
MInt d = 0, l = 0; d < noDomains(); d++) {
2368 MPI_Group tmp_group;
2371 MPI_Group_incl(tmp_group, noRanks, ownerList.getPointer(), &group, AT_);
2390 name <<
"out/" << domainId() <<
"_boundary.vtk";
2392 st.open(name.str());
2393 st <<
"# vtk DataFile Version 3.0" << endl;
2394 st <<
"vtk output" << endl;
2395 st <<
"ASCII" << endl;
2396 st <<
"DATASET POLYDATA" << endl;
2397 st <<
"POINTS " << num <<
" float" << endl;
2399 for(
MInt i = 0; i < num; i++) {
2400 for(
MInt d = 0; d < 3; d++)
2401 st << bndVs[i][d] <<
" ";
2405 st <<
"LINES " << 1 <<
" " << num + 2 << endl;
2406 st << num + 1 <<
" ";
2407 for(
MInt i = 0; i < num; i++) {
2425 constexpr MFloat L_ref = 1.0;
2426 const MFloat dx = c_cellLengthAtLevel(maxLevel());
2427 m_dt = dx * m_Ma * LBCS / L_ref;
2428 m_time = getCurrentTimeStep() * m_dt;
2448 m_log <<
" * segment owned by: ";
2449 m_geometry->determineSegmentOwnership(segmentId, &own, &sumowners, &firstOwner, owners.getPointer());
2450 for(
MInt d = 0; d < noDomains(); d++)
2451 if(owners[d] > 0)
m_log << d <<
" ";
2454 m_log <<
" * sum of owners: " << sumowners << endl;
2455 m_log <<
" * root of communication: " << firstOwner << endl;
2459 if(sumowners == 1) {
2461 if(own) result = calcCharLenAll(segmentId);
2463 MPI_Bcast(&result, 1, MPI_DOUBLE, firstOwner, mpiComm(), AT_,
"result");
2474 createMPIComm(owners.getPointer(), sumowners, &charComm);
2478 MInt offStart = m_geometry->m_segmentOffsets[segmentId];
2479 MInt offEnd = m_geometry->m_segmentOffsets[segmentId + 1];
2480 MInt numElements = offEnd - offStart;
2482 m_log <<
" * number of local triangles: " << numElements << endl;
2483 m_log <<
" * segment offsets: " << m_geometry->m_segmentOffsets[segmentId] <<
" - "
2484 << m_geometry->m_segmentOffsets[segmentId + 1] << endl;
2487 MInt noTriInfo = nDim * nDim;
2491 for(
MInt t = offStart, i = 0, j = 0; t < offEnd; t++, i++) {
2492 myOriginalIds[i] = m_geometry->elements[t].m_originalId;
2493 for(
MInt v = 0; v < nDim; v++)
2494 for(
MInt d = 0; d < nDim; d++, j++)
2495 segTriangles[j] = m_geometry->elements[t].m_vertices[v][d];
2499 MPI_Allgather(&numElements, 1, MPI_INT, numElemPerCPU.getPointer(), 1, MPI_INT, charComm, AT_,
"numElements",
2500 "numElemPerCPU.getPointer()");
2503 m_log <<
" * triangles per domain: ";
2504 MInt sumallelem = 0;
2505 for(
MInt i = 0; i < sumowners; i++) {
2506 m_log << numElemPerCPU[i] <<
" ";
2507 sumallelem += numElemPerCPU[i];
2508 numTriInfoPerCPU[i] = numElemPerCPU[i] * noTriInfo;
2511 m_log <<
" * sum of global triangles: " << sumallelem << endl;
2518 for(
MInt d = 1; d < sumowners; d++) {
2519 displOrig[d] = displOrig[d - 1] + numElemPerCPU[d - 1];
2520 displTris[d] = displTris[d - 1] + numTriInfoPerCPU[d - 1];
2524 MPI_Gatherv(myOriginalIds.getPointer(), numElements, MPI_INT, allOriginalIds.getPointer(),
2525 numElemPerCPU.getPointer(), displOrig.getPointer(), MPI_INT, 0, charComm, AT_,
2526 "myOriginalIds.getPointer()",
"allOriginalIds.getPointer()");
2529 MPI_Gatherv(segTriangles.getPointer(), noTriInfo * numElements, MPI_DOUBLE, allSegTriangles.getPointer(),
2530 numTriInfoPerCPU.getPointer(), displTris.getPointer(), MPI_DOUBLE, 0, charComm, AT_,
2531 "segTriangles.getPointer()",
"allSegTriangles.getPointer()");
2534 if(domainId() == firstOwner) {
2535 set<MInt> uniqueTriangles;
2536 for(
MInt i = 0; i < sumallelem; i++)
2537 uniqueTriangles.insert(allOriginalIds[i]);
2539 MInt noUniqueTris = uniqueTriangles.size();
2540 m_log <<
" * sum of unique triangles: " << noUniqueTris << endl;
2543 for(
MInt i = 0; i < noUniqueTris; i++)
2548 for(
MInt i = 0, j = 0; i < sumallelem; i++) {
2549 MInt dist =
distance(uniqueTriangles.begin(), uniqueTriangles.find(allOriginalIds[i]));
2550 if(
dist != noUniqueTris && dbl[
dist] == 0) {
2551 keepOffsets[j] = i * noTriInfo;
2558 MFloat** bndVs = m_geometry->GetBoundaryVertices(segmentId, allSegTriangles.getPointer(),
2559 keepOffsets.getPointer(), noUniqueTris, &num);
2561 if(m_geometry->m_debugParGeom) writeSegmentBoundaryVTK(bndVs, num);
2564 MFloat circ = m_geometry->calcCircumference(bndVs, num);
2565 m_log <<
" * circumference: " << circ << endl;
2568 MFloat size = m_geometry->GetBoundarySize(allSegTriangles.getPointer(), keepOffsets.getPointer(), noUniqueTris);
2569 m_log <<
" * area: " << size << endl;
2572 MFloat hydraulic_diam = 4 * size / circ;
2573 m_log <<
" * hydraulic diameter: " << hydraulic_diam << endl;
2576 std::array<MFloat, nDim * 2> bBox;
2577 m_geometry->getBoundingBox(bBox.data());
2582 for(
MInt i = 0; i < nDim; i++)
2583 if(fabs(bBox[i + nDim] - bBox[i]) > maxlength) maxlength = fabs(bBox[i + nDim] - bBox[i]);
2585 diameter = hydraulic_diam / (maxlength /
FPOW2(maxLevel())) / reductionFactor();
2588 MPI_Bcast(&diameter, 1, MPI_DOUBLE, firstOwner, mpiComm(), AT_,
"diameter");
2622 MFloat** bndVs = m_geometry->GetBoundaryVertices(segmentId,
nullptr,
nullptr, 0, &num);
2625 MFloat circ = m_geometry->calcCircumference(bndVs, num);
2628 MFloat size = m_geometry->GetBoundarySize(segmentId);
2631 MFloat hydraulic_diam = 4 * size / circ;
2634 std::array<MFloat, nDim * 2> bBox;
2635 m_geometry->getBoundingBox(bBox.data());
2640 for(
MInt i = 0; i < nDim; i++)
2641 if(fabs(bBox[i + nDim] - bBox[i]) > maxlength) maxlength = fabs(bBox[i + nDim] - bBox[i]);
2643 m_log <<
" * circumference: " << circ << endl;
2644 m_log <<
" * area: " << size << endl;
2645 m_log <<
" * hydraulic diameter: " << hydraulic_diam << endl;
2647 return hydraulic_diam / (maxlength /
FPOW2(maxLevel())) / reductionFactor();
2662 m_log <<
" + Calculating characteristic length:" << endl;
2663 m_log <<
" - type: " << (m_geometry->m_parallelGeometry ?
"parallel" :
"serial") << endl;
2664 m_log <<
" - segment id: " << segmentId << endl;
2666 if constexpr(nDim == 2) {
2667 stringstream errorMsg;
2668 errorMsg <<
"ERROR: no implementation for the calulation of the characteristic length in 2D!" << endl;
2669 m_log << errorMsg.str();
2670 TERMM(1, errorMsg.str());
2674 if(m_geometry->m_parallelGeometry)
2675 ret = calcCharLenParallelSplit(segmentId);
2677 ret = calcCharLenAll(segmentId);
2679 m_log <<
" * charcteristic length: " << ret <<
" [cell units]" << endl;
2688 static ofstream ofl;
2693 ofl.open(
"cellInfo.dat", ios_base::app);
2695 ofl <<
" cell informations: " << endl;
2696 ofl <<
" -------------------" << endl;
2697 ofl <<
" id=" <<
cellId <<
", domain=" << domainId() <<
", level=" << a_level(cellId) << endl;
2698 ofl <<
" (x,y,z)=(" << a_coordinate(cellId, 0) <<
", " << a_coordinate(cellId, 1) <<
", " << a_coordinate(cellId, 2)
2700 ofl <<
" neighbors: " << endl;
2701 for(
MInt k = 0; k < 26; k++) {
2702 nghbrId = c_neighborId(cellId, k);
2703 ofl <<
" " << k <<
": noNghbrs=" << a_hasNeighbor(cellId, k) <<
", id=" << nghbrId << endl;
2704 ofl <<
" (x,y,z)=(" << a_coordinate(nghbrId, 0) <<
", " << a_coordinate(nghbrId, 1) <<
", "
2705 << a_coordinate(nghbrId, 2) <<
")" << endl;
2707 ofl <<
" children: " << endl;
2708 ofl <<
" noChildren: " << c_noChildren(cellId) << endl;
2710 ofl <<
" " << k <<
": " << c_childId(cellId, k) << endl;
2712 ofl <<
" parent: " << c_parentId(cellId) << endl;
2714 if(a_isInterface(cellId)) {
2715 for(
MInt k = 0; k < 26; k++) {
2716 nghbrId = c_neighborId(c_parentId(cellId), k);
2717 ofl <<
" " << k <<
": parent noNghbrs=" << a_hasNeighbor(c_parentId(cellId), k) <<
", id=" << nghbrId
2719 ofl <<
" (x,y,z)=(" << a_coordinate(nghbrId, 0) <<
", " << a_coordinate(nghbrId, 1) <<
", "
2720 << a_coordinate(nghbrId, 2) <<
")" << endl;
2732 MInt tmpCounter = 0;
2738 if constexpr(nDim == 3) {
2740 ofl <<
" TITLE = \" AMAZONAS Tecfile \" " << endl;
2741 ofl << R
"( VARIABLES = "X", "Y", "Z", "U", "V", "W", "RHO")" << endl;
2742 ofl << " ZONE N=" << m_gridPoints->size() <<
", ";
2743 ofl <<
" E=" << m_extractedCells->size() <<
", ";
2747 ofl <<
" DATAPACKING="
2748 <<
"SOLVER" << endl;
2749 ofl <<
" VARLOCATION=([4, 5, 6, 7]=CELLCENTERED)" << endl;
2753 ofl <<
" TITLE = \" AMAZONAS Tecfile \" " << endl;
2754 ofl << R
"( VARIABLES = "X", "Y", "U", "V", "RHO")" << endl;
2755 ofl << " ZONE N=" << m_gridPoints->size() <<
", ";
2756 ofl <<
" E=" << m_extractedCells->size() <<
", ";
2758 <<
"FEQUADRILATERAL"
2760 ofl <<
" DATAPACKING="
2761 <<
"SOLVER" << endl;
2762 ofl <<
" VARLOCATION=([3, 4, 5]=CELLCENTERED)" << endl;
2766 for(
MInt axisId = 0; axisId < nDim; axisId++) {
2767 for(
MInt pointId = 0; pointId < m_gridPoints->size(); pointId++) {
2768 ofl << m_gridPoints->a[pointId].m_coordinates[axisId] <<
" ";
2770 if(tmpCounter == 10) {
2779 for(
MInt varId = 0; varId < nDim + 1; varId++) {
2781 ofl << a_variable(m_extractedCells->a[cellId].m_cellId, varId) <<
" ";
2784 if(tmpCounter == 8) {
2799 ofl << m_extractedCells->a[
cellId].m_pointIds[2] + 1 <<
" ";
2800 ofl << m_extractedCells->a[
cellId].m_pointIds[3] + 1 <<
" ";
2801 ofl << m_extractedCells->a[
cellId].m_pointIds[1] + 1 <<
" ";
2802 ofl << m_extractedCells->a[
cellId].m_pointIds[0] + 1 <<
" ";
2803 if constexpr(nDim == 3) {
2804 ofl << m_extractedCells->a[
cellId].m_pointIds[6] + 1 <<
" ";
2805 ofl << m_extractedCells->a[
cellId].m_pointIds[7] + 1 <<
" ";
2806 ofl << m_extractedCells->a[
cellId].m_pointIds[5] + 1 <<
" ";
2807 ofl << m_extractedCells->a[
cellId].m_pointIds[4] + 1 <<
" ";
2813 cerr <<
"ERROR! COULD NOT OPEN FILE " << fileName <<
" for writing! " << endl;
2823 MInt noFieldElements = nDim + 1;
2825#define WRITE_CELLIDS
2826#define WRITE_OLD_VARIABLES
2827#define WRITE_DISTRIBUTIONS
2831#ifdef WRITE_OLD_VARIABLES
2832 noFieldElements += nDim + 1;
2834#ifdef WRITE_DISTRIBUTIONS
2835 noFieldElements += m_noDistributions;
2843 ofl <<
"# vtk DataFile Version 3.0 " << endl;
2844 ofl <<
"Amazonas output file" << endl;
2845 ofl <<
"BINARY " << endl;
2846 ofl <<
"DATASET UNSTRUCTURED_GRID " << endl;
2847 ofl <<
"POINTS " << m_gridPoints->size() <<
" float" << endl;
2851 float tmpCoordinate = 0;
2852 for(
MInt pointId = 0; pointId < m_gridPoints->size(); pointId++) {
2853 for(
MInt axisId = 0; axisId < nDim; axisId++) {
2854 tmpCoordinate = (float)m_gridPoints->a[pointId].m_coordinates[axisId];
2856 tmpCoordinate =
floatSwap(tmpCoordinate);
2858 ofl.write(
reinterpret_cast<char*
>(&tmpCoordinate),
sizeof(float));
2863 if constexpr(nDim == 2) {
2867 ofl.write(
reinterpret_cast<char*
>(&tmpCoordinate),
sizeof(
float));
2871 ofl.open(fileName, ios_base::app);
2875 ofl <<
"CELLS " << m_extractedCells->size() <<
" " << m_extractedCells->size() * (number + 1) << endl;
2879 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
2886 ofl.write(
reinterpret_cast<char*
>(&number),
sizeof(
MInt));
2888 for(
MInt points = 0; points <
IPOW2(nDim); points++) {
2889 pointId = m_extractedCells->a[
cellId].m_pointIds[points];
2893 ofl.write(
reinterpret_cast<char*
>(&pointId),
sizeof(
MInt));
2897 ofl.open(fileName, ios_base::app);
2899 ofl <<
"CELL_TYPES " << m_extractedCells->size() << endl;
2901 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
2902 if constexpr(nDim == 3) {
2912 ofl.write(
reinterpret_cast<char*
>(&pointId),
sizeof(
MInt));
2916 ofl.open(fileName, ios_base::app);
2918 ofl <<
"CELL_DATA " << m_extractedCells->size() << endl;
2919 ofl <<
"SCALARS scalars double" << endl;
2920 ofl <<
"LOOKUP_TABLE default" << endl;
2923 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
2929 ofl.write(
reinterpret_cast<char*
>(&fnumber),
sizeof(
MFloat));
2933#ifdef WRITE_VECTOR_DATA
2935 ofl.open(fileName, ios_base::app);
2937 ofl <<
"VECTORS vectors double" << endl;
2940 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
2942 for(
MInt varId = 0; varId < nDim; varId++) {
2943 fnumber = a_variable(m_extractedCells->a[cellId].m_cellId, varId);
2947 ofl.write(
reinterpret_cast<char*
>(&fnumber),
sizeof(
MFloat));
2955 ofl.open(fileName, ios_base::app);
2957 ofl <<
"FIELD FieldData " << noFieldElements << endl;
2959 for(
MInt varId = 0; varId < nDim + 1; varId++) {
2961 ofl.open(fileName, ios_base::app);
2963 ofl <<
"variables" << varId <<
" 1 " << m_extractedCells->size() <<
" double" << endl;
2965 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
2973 fnumber = a_variable(m_extractedCells->a[cellId].m_cellId, varId);
2977 ofl.write(
reinterpret_cast<char*
>(&fnumber),
sizeof(
MFloat));
2981#ifdef WRITE_OLD_VARIABLES
2982 for(
MInt varId = 0; varId < nDim + 1; varId++) {
2984 ofl.open(fileName, ios_base::app);
2986 ofl <<
"residual" << varId <<
" 1 " << m_extractedCells->size() <<
" double" << endl;
2988 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
2990 fnumber = a_variable(m_extractedCells->a[cellId].m_cellId, varId)
2991 - a_oldVariable(m_extractedCells->a[cellId].m_cellId, varId);
2995 ofl.write(
reinterpret_cast<char*
>(&fnumber),
sizeof(
MFloat));
2999#ifdef WRITE_DISTRIBUTIONS
3001 for(
MInt distId = 0; distId < m_noDistributions; distId++) {
3003 ofl.open(fileName, ios_base::app);
3005 ofl <<
"distributions" << distId <<
" 1 " << m_extractedCells->size() <<
" double" << endl;
3007 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
3009 fnumber = a_distribution(m_extractedCells->a[cellId].m_cellId, distId);
3013 ofl.write(
reinterpret_cast<char*
>(&fnumber),
sizeof(
MFloat));
3022 ofl.open(fileName, ios_base::app);
3025 <<
" 1 " << m_extractedCells->size() <<
" int" << endl;
3027 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
3029 inumber = m_extractedCells->a[
cellId].m_cellId;
3033 ofl.write(
reinterpret_cast<char*
>(&inumber),
sizeof(
MInt));
3039 cerr <<
"ERROR! COULD NOT OPEN FILE " << fileName <<
" for writing! " << endl;
3047 MInt tmpCounter = 0;
3052 if constexpr(nDim == 3) {
3054 ofl <<
"# vtk DataFile Version 2.0 " << endl;
3055 ofl <<
"Amazonas output file" << endl;
3056 ofl <<
"ASCII " << endl;
3057 ofl <<
"DATASET UNSTRUCTURED_GRID " << endl;
3058 ofl <<
"POINTS " << m_gridPoints->size() <<
" float" << endl;
3060 for(
MInt pointId = 0; pointId < m_gridPoints->size(); pointId++) {
3061 for(
MInt axisId = 0; axisId < nDim; axisId++) {
3062 ofl << m_gridPoints->a[pointId].m_coordinates[axisId] <<
" ";
3066 if(tmpCounter == 10) {
3074 ofl <<
"CELLS " << m_extractedCells->size() <<
" " << m_extractedCells->size() * 9 << endl;
3079 ofl << m_extractedCells->a[
cellId].m_pointIds[0] <<
" ";
3080 ofl << m_extractedCells->a[
cellId].m_pointIds[1] <<
" ";
3081 ofl << m_extractedCells->a[
cellId].m_pointIds[2] <<
" ";
3082 ofl << m_extractedCells->a[
cellId].m_pointIds[3] <<
" ";
3083 if constexpr(nDim == 3) {
3084 ofl << m_extractedCells->a[
cellId].m_pointIds[4] <<
" ";
3085 ofl << m_extractedCells->a[
cellId].m_pointIds[5] <<
" ";
3086 ofl << m_extractedCells->a[
cellId].m_pointIds[6] <<
" ";
3087 ofl << m_extractedCells->a[
cellId].m_pointIds[7] <<
" ";
3092 ofl <<
"CELL_TYPES " << m_extractedCells->size() << endl;
3094 ofl <<
"11 " << endl;
3096 ofl <<
"CELL_DATA " << m_extractedCells->size() << endl;
3097 ofl <<
"SCALARS scalars float" << endl;
3098 ofl <<
"LOOKUP_TABLE default" << endl;
3100 ofl <<
"1.0" << endl;
3102 ofl <<
"FIELD FieldData " << nDim + 1 << endl;
3104 for(
MInt varId = 0; varId < nDim + 1; varId++) {
3105 ofl <<
"variables" << varId <<
" 1 " << m_extractedCells->size() <<
" float" << endl;
3108 ofl << a_variable(m_extractedCells->a[cellId].m_cellId, varId) <<
" ";
3116 ofl <<
"FIELD FieldData " << nDim + 1 << endl;
3117 for(
MInt varId = 0; varId < nDim + 1; varId++) {
3118 ofl <<
"oldVariables" << varId <<
" 1 " << m_extractedCells->size() <<
" float" << endl;
3121 ofl << a_oldVariable(m_extractedCells->a[cellId].m_cellId, varId) <<
" ";
3130 cerr <<
"ERROR! COULD NOT OPEN FILE " << fileName <<
" for writing! " << endl;
3159 if constexpr(nDim != 2 && nDim != 3) {
3160 cerr <<
" In global function loadGridFlowVariablesThermalPar: wrong number of dimensions !" << endl;
3166 ParallelIo parallelIo(fileName, PIO_READ, mpiComm());
3169 if(!m_initFromRestartFile) {
3170 parallelIo.getAttribute(&m_time,
"time");
3174 if(m_velocityControl.restart) {
3175 parallelIo.getAttribute(&m_velocityControl.lastGlobalAvgV,
"lastGlobalAvgV");
3176 parallelIo.getAttribute(&m_velocityControl.previousError,
"previousError");
3177 parallelIo.getAttribute(&m_velocityControl.integratedError,
"integratedError");
3178 parallelIo.getAttribute(&m_velocityControl.derivedError,
"derivedError");
3179 for(
MInt i = 0; i < nDim; i++) {
3180 parallelIo.getAttribute(&m_volumeAccel[i],
"volumeAcceleration_" + std::to_string(i));
3184 MInt coarse_count = 0;
3185 parallelIo.readScalar(&coarse_count,
"noCells");
3188 MInt local_count = 0;
3189 for(
MInt i = 0; i < noInternalCells(); i++) {
3190 if(c_noChildren(i) > 0) local_count++;
3194 for(
MInt d = 0; d < noDomains(); d++)
3196 nonFine_count[d] = local_count;
3198 nonFine_count[d] = 0;
3200 MPI_Allreduce(MPI_IN_PLACE, nonFine_count.getPointer(), noDomains(), MPI_INT, MPI_SUM, mpiComm(), AT_,
"MPI_IN_PLACE",
3201 "nonFine_count.getPointer()");
3204 for(
MInt d = 0; d < domainId(); d++)
3205 off += nonFine_count[d];
3208 for(
MInt d = 0; d < noDomains(); d++)
3209 all += nonFine_count[d];
3211 m_log <<
" - loading information: " << endl;
3212 m_log <<
" * number of all cells: " << all << endl;
3213 m_log <<
" * number of local cells: " << local_count << endl;
3214 m_log <<
" * domain offset: " << off << endl;
3216 const MPI_Offset firstGlobalId = off;
3217 const MPI_Offset localNoCells = local_count;
3218 parallelIo.setOffset(localNoCells, firstGlobalId);
3220 m_log <<
" - determining ids of the coarse cells" << endl;
3223 for(
MInt i = 0, j = 0; i < noInternalCells(); ++i)
3224 if(c_noChildren(i) > 0) {
3229 m_log <<
" - reading variables " << endl;
3235 for(
MInt v = 0; v < m_noVariables; v++) {
3236 MString name =
"variables" + to_string(v);
3238 parallelIo.readArray(tmparray.getPointer(), name);
3239 for(
MInt i = 0; i < local_count; ++i)
3240 a_variable(coarse_ids[i], v) = tmparray.p[i];
3244 for(
MInt v = 0; v < m_noVariables; v++) {
3245 MString name =
"oldVariables" + to_string(v);
3247 parallelIo.readArray(tmparray.getPointer(), name);
3248 for(
MInt i = 0; i < local_count; ++i)
3249 a_variable(coarse_ids[i], v) = tmparray.p[i];
3253 m_log <<
" - exchanging coarse window cells" << endl;
3256 m_log <<
" * collect windows and halos on a coarse level " << endl;
3257 vector<vector<MInt>> coarseWindowPD;
3258 vector<vector<MInt>> coarseHaloPD;
3259 for(
MInt d = 0; d < noNeighborDomains(); d++) {
3261 for(
MInt w = 0; w < noWindowCells(d); w++)
3262 if(a_level(windowCell(d, w)) == maxLevel() - 1) cW.push_back(windowCell(d, w));
3265 for(
MInt h = 0; h < noHaloCells(d); h++)
3266 if(a_level(haloCell(d, h)) == maxLevel() - 1) cH.push_back(haloCell(d, h));
3268 coarseWindowPD.push_back(cW);
3269 coarseHaloPD.push_back(cH);
3272 m_log <<
" = number of coarse window cells for " << noNeighborDomains() <<
" neighboring domains" << endl;
3273 for(
MInt d = 0; d < noNeighborDomains(); d++)
3274 m_log <<
" # domain " << neighborDomain(d) <<
":" << coarseWindowPD[d].size() << endl;
3275 m_log <<
" = number of coarse halo cells for " << noNeighborDomains() <<
" neighboring domains" << endl;
3276 for(
MInt d = 0; d < noNeighborDomains(); d++)
3277 m_log <<
" # domain " << neighborDomain(d) <<
":" << coarseHaloPD[d].size() << endl;
3280 m_log <<
" * constructing send buffer with size: ";
3281 MInt totalWindowBuf = 0;
3283 for(
MInt d = 0; d < noNeighborDomains(); d++) {
3284 windowOff[d] = totalWindowBuf;
3285 totalWindowBuf += coarseWindowPD[d].size() * 2 * m_noVariables;
3288 m_log << totalWindowBuf << endl;
3291 for(
MInt d = 0; d < noNeighborDomains(); d++) {
3292 MInt pos = windowOff[d];
3293 for(
MInt c = 0; c < (
MInt)coarseWindowPD[d].size(); c++) {
3294 for(
MInt v = 0; v < m_noVariables; v++, pos++)
3295 sendBuf[pos] = a_variable(coarseWindowPD[d][c], v);
3296 for(
MInt v = 0; v < m_noVariables; v++, pos++)
3297 sendBuf[pos] = a_oldVariable(coarseWindowPD[d][c], v);
3299 MPI_Issend(&sendBuf[windowOff[d]], coarseWindowPD[d].size() * 2 * m_noVariables, MPI_DOUBLE, neighborDomain(d), 0,
3300 mpiComm(), &mpi_request[d], AT_,
"sendBuf[windowOff[d]]");
3304 m_log <<
" * constructing receive buffer and receive" << endl;
3305 MInt totalHaloBuf = 0;
3307 for(
MInt d = 0; d < noNeighborDomains(); d++) {
3308 haloOff[d] = totalHaloBuf;
3309 totalHaloBuf += coarseHaloPD[d].size() * 2 * m_noVariables;
3315 for(
MInt d = 0; d < noNeighborDomains(); d++)
3316 MPI_Recv(&receiveBuf[haloOff[d]], coarseHaloPD[d].size() * 2 * m_noVariables, MPI_DOUBLE, neighborDomain(d), 0,
3317 mpiComm(), &status, AT_,
"receiveBuf[haloOff[d]]");
3319 for(
MInt d = 0; d < noNeighborDomains(); d++)
3320 MPI_Wait(&mpi_request[d], &status, AT_);
3323 m_log <<
" * distributing received information" << endl;
3324 for(
MInt d = 0; d < noNeighborDomains(); d++) {
3325 MInt pos = haloOff[d];
3326 for(
MInt c = 0; c < (
MInt)coarseHaloPD[d].size(); c++) {
3327 for(
MInt v = 0; v < m_noVariables; v++, pos++)
3328 a_variable(coarseHaloPD[d][c], v) = receiveBuf[pos];
3329 for(
MInt v = 0; v < m_noVariables; v++, pos++)
3330 a_oldVariable(coarseHaloPD[d][c], v) = receiveBuf[pos];
3335 m_log <<
" - filling all finer cells" << endl;
3339 if constexpr(nDim == 2) {
3340 for(
MInt d = 0; d < elems; d++) {
3341 dirs[d] = childPos2D[d];
3344 for(
MInt d = 0; d < elems; d++) {
3345 dirs[d] = childPos3D[d];
3350 m_log <<
" * filling all cells with nieghbors" << endl;
3351 for(
MInt i = 0; i < noInternalCells(); i++)
3352 if(a_level(i) == maxLevel()) {
3353 MInt parent = c_parentId(i);
3357 for(
MInt c = 0; c < elems; c++)
3358 if(c_childId(parent, c) == i) {
3366 if(a_hasNeighbor(parent, dir)) {
3367 MInt parentNeighbor = c_neighborId(parent, dir);
3368 for(
MInt v = 0; v < m_noVariables; v++) {
3369 a_variable(i, v) = F3B4 * a_variable(parent, v) + F1B4 * a_variable(parentNeighbor, v);
3370 a_oldVariable(i, v) = F3B4 * a_oldVariable(parent, v) + F1B4 * a_oldVariable(parentNeighbor, v);
3375 for(
MInt v = 0; v < m_noVariables; v++) {
3376 a_variable(i, v) = a_variable(parent, v);
3377 a_oldVariable(i, v) = a_oldVariable(parent, v);
3383 m_log <<
" * filling halos" << endl;
3385 MInt noNeigh = pow(nDim, nDim);
3386 for(
MInt d = 0; d < noNeighborDomains(); d++)
3387 for(
MInt h = 0; h < noHaloCells(d); h++) {
3388 MInt haloId = haloCell(d, h);
3390 if(a_level(haloId) == maxLevel()) {
3392 MInt haloWinNeigh = -1;
3394 for(; dir < noNeigh; dir++)
3395 if(a_hasNeighbor(haloId, dir) && a_isWindow(c_neighborId(haloId, dir))) {
3396 haloWinNeigh = c_neighborId(haloId, dir);
3401 if(haloWinNeigh >= 0) {
3402 MInt winParent = c_parentId(haloWinNeigh);
3408 if(winParent >= 0 && a_hasNeighbor(winParent, opp)) parent = c_neighborId(winParent, opp);
3410 if(parent < 0) parent = winParent;
3412 for(
MInt v = 0; v < m_noVariables; v++) {
3413 a_variable(haloId, v) = a_variable(parent, v);
3414 a_oldVariable(haloId, v) = a_oldVariable(parent, v);
3436 if constexpr(nDim != 2 && nDim != 3) {
3437 cerr <<
" In global function loadGridFlowVariablesThermalPar: wrong number of dimensions !" << endl;
3443 ParallelIo parallelIo(fileName, PIO_READ, mpiComm());
3446 if(!m_initFromRestartFile) {
3447 parallelIo.getAttribute(&m_time,
"time");
3451 if(m_velocityControl.restart) {
3452 parallelIo.getAttribute(&m_velocityControl.lastGlobalAvgV,
"lastGlobalAvgV");
3453 parallelIo.getAttribute(&m_velocityControl.previousError,
"previousError");
3454 parallelIo.getAttribute(&m_velocityControl.integratedError,
"integratedError");
3455 parallelIo.getAttribute(&m_velocityControl.derivedError,
"derivedError");
3456 for(
MInt i = 0; i < nDim; i++) {
3457 parallelIo.getAttribute(&m_volumeAccel[i],
"volumeAcceleration_" + std::to_string(i));
3461 MInt pv_veloc[nDim];
3462 if constexpr(nDim == 2) {
3463 pv_veloc[0] = PV->U;
3464 pv_veloc[1] = PV->V;
3466 pv_veloc[0] = PV->U;
3467 pv_veloc[1] = PV->V;
3468 pv_veloc[2] = PV->W;
3470 MInt pvrho = PV->RHO;
3483 MString oldDistributionsThermal;
3484 MString distributionsTransport;
3485 MString oldDistributionsTransport;
3486 if constexpr(nDim == 3) {
3487 rhoName =
"variables3";
3488 oldrhoName =
"oldVariables3";
3489 temperatureName =
"variables4";
3490 oldtemperatureName =
"oldVariables4";
3491 concentrationName =
"variables5";
3492 oldconcentrationName =
"oldVariables5";
3493 EELiquidName =
"oldVariables4";
3495 rhoName =
"variables2";
3496 oldrhoName =
"oldVariables2";
3497 temperatureName =
"variables3";
3498 oldtemperatureName =
"oldVariables3";
3499 concentrationName =
"variables4";
3500 oldconcentrationName =
"oldVariables4";
3505 const MPI_Offset firstGlobalId = domainOffset(domainId());
3506 const MPI_Offset localNoCells = noInternalCells();
3507 parallelIo.setOffset(localNoCells, firstGlobalId);
3515 parallelIo.readArray(tmparray.getPointer(),
"variables0");
3516 for(
MInt i = 0; i < noInternalCells(); ++i)
3517 a_variable(i, pv_veloc[0]) = tmparray[i];
3520 parallelIo.readArray(tmparray.getPointer(),
"oldVariables0");
3521 for(
MInt i = 0; i < noInternalCells(); ++i)
3522 a_oldVariable(i, pv_veloc[0]) = tmparray[i];
3525 parallelIo.readArray(tmparray.getPointer(),
"variables1");
3526 for(
MInt i = 0; i < noInternalCells(); ++i)
3527 a_variable(i, pv_veloc[1]) = tmparray.p[i];
3530 parallelIo.readArray(tmparray.getPointer(),
"oldVariables1");
3531 for(
MInt i = 0; i < noInternalCells(); ++i)
3532 a_oldVariable(i, pv_veloc[1]) = tmparray.p[i];
3535 parallelIo.readArray(tmparray.getPointer(), rhoName);
3536 for(
MInt i = 0; i < noInternalCells(); ++i)
3537 a_variable(i, pvrho) = tmparray.p[i];
3540 parallelIo.readArray(tmparray.getPointer(), oldrhoName);
3541 for(
MInt i = 0; i < noInternalCells(); ++i)
3542 a_oldVariable(i, pvrho) = tmparray.p[i];
3546 parallelIo.readArray(tmparray.getPointer(), temperatureName);
3547 for(
MInt i = 0; i < noInternalCells(); ++i)
3548 a_variable(i, pvt) = tmparray.p[i];
3551 parallelIo.readArray(tmparray.getPointer(), oldtemperatureName);
3552 for(
MInt i = 0; i < noInternalCells(); ++i)
3553 a_oldVariable(i, pvt) = tmparray.p[i];
3557 parallelIo.readArray(tmparray.getPointer(), concentrationName);
3558 for(
MInt i = 0; i < noInternalCells(); ++i)
3559 a_variable(i, pvc) = tmparray.p[i];
3562 parallelIo.readArray(tmparray.getPointer(), oldconcentrationName);
3563 for(
MInt i = 0; i < noInternalCells(); ++i)
3564 a_oldVariable(i, pvc) = tmparray.p[i];
3567 if(!m_EELiquid.restartWithoutAlpha) {
3568 parallelIo.readArray(tmparray.getPointer(), EELiquidName);
3569 for(
MInt i = 0; i < noInternalCells(); ++i)
3570 a_alphaGas(i) = tmparray.p[i];
3572 for(
MInt i = 0; i < noInternalCells(); ++i)
3573 a_alphaGas(i) = 0.0;
3579 if constexpr(nDim == 3) {
3583 parallelIo.readArray(tmparray.getPointer(),
"variables2");
3584 for(
MInt i = 0; i < noInternalCells(); ++i)
3585 a_variable(i, pv_veloc[2]) = tmparray.p[i];
3588 parallelIo.readArray(tmparray.getPointer(),
"oldVariables2");
3589 for(
MInt i = 0; i < noInternalCells(); ++i)
3590 a_oldVariable(i, pv_veloc[2]) = tmparray.p[i];
3595 for(
MInt j = 0; j < m_noDistributions; j++) {
3596 distributions =
"distributions" + to_string(j);
3599 parallelIo.readArray(tmp_distributions.getPointer(), distributions);
3600 for(
MInt i = 0; i < noInternalCells(); ++i)
3601 a_distribution(i, j) = tmp_distributions.p[i];
3605 for(
MInt j = 0; j < m_noDistributions; j++) {
3606 oldDistributions =
"oldDistributions" + to_string(j);
3608 MFloatScratchSpace tmp_oldDistributions(noInternalCells(), AT_,
"tmp_oldDistributions");
3609 parallelIo.readArray(tmp_oldDistributions.getPointer(), oldDistributions);
3610 for(
MInt i = 0; i < noInternalCells(); ++i)
3611 a_oldDistribution(i, j) = tmp_oldDistributions.p[i];
3616 for(
MInt j = 0; j < m_noDistributions; j++) {
3617 distributionsThermal =
"distributionsThermal" + to_string(j);
3619 MFloatScratchSpace tmp_distributionsThermal(noInternalCells(), AT_,
"tmp_distributionsThermal");
3620 parallelIo.readArray(tmp_distributionsThermal.getPointer(), distributionsThermal);
3621 for(
MInt i = 0; i < noInternalCells(); ++i)
3622 a_distributionThermal(i, j) = tmp_distributionsThermal.p[i];
3626 for(
MInt j = 0; j < m_noDistributions; j++) {
3627 oldDistributionsThermal =
"oldDistributionsThermal" + to_string(j);
3629 MFloatScratchSpace tmp_oldDistributionsThermal(noInternalCells(), AT_,
"tmp_oldDistributionsThermal");
3630 parallelIo.readArray(tmp_oldDistributionsThermal.getPointer(), oldDistributionsThermal);
3631 for(
MInt i = 0; i < noInternalCells(); ++i)
3632 a_oldDistributionThermal(i, j) = tmp_oldDistributionsThermal.p[i];
3635 for(
MInt i = 0; i < a_noCells(); i++) {
3636 a_kappa(i) = m_kappa;
3642 for(
MInt j = 0; j < m_noDistributions; j++) {
3643 distributionsTransport =
"distributionsTransport" + to_string(j);
3645 MFloatScratchSpace tmp_distributionsTransport(noInternalCells(), AT_,
"tmp_distributionsTransport");
3646 parallelIo.readArray(tmp_distributionsTransport.getPointer(), distributionsTransport);
3647 for(
MInt i = 0; i < noInternalCells(); ++i)
3648 a_distributionTransport(i, j) = tmp_distributionsTransport.p[i];
3652 for(
MInt j = 0; j < m_noDistributions; j++) {
3653 oldDistributionsTransport =
"oldDistributionsTransport" + to_string(j);
3655 MFloatScratchSpace tmp_oldDistributionsTransport(noInternalCells(), AT_,
"tmp_oldDistributionsTransport");
3656 parallelIo.readArray(tmp_oldDistributionsTransport.getPointer(), oldDistributionsTransport);
3657 for(
MInt i = 0; i < noInternalCells(); ++i)
3658 a_oldDistributionTransport(i, j) = tmp_oldDistributionsTransport.p[i];
3661 for(
MInt i = 0; i < a_noCells(); i++) {
3662 a_diffusivity(i) = m_diffusivity;
3666 if(!m_nonNewtonian) {
3667 for(
MInt i = 0; i < a_noCells(); i++) {
3672 parallelIo.readArray(tmp_nu.getPointer(),
"viscosity");
3673 for(
MInt i = 0; i < noInternalCells(); i++) {
3674 initNu(i, tmp_nu(i));
3679 m_bndCnd->bcDataReadRestartData(parallelIo);
3682 if(noDomains() > 1) {
3684 exchangeOldDistributions();
3685 }
else if(noNeighborDomains() > 0) {
3688 std::stringstream ss;
3689 ss <<
"WARNING: Serial run and periodic (?) is problematic as halo cells are"
3690 " not exchanged even though periodic BC make usage of it!"
3693 std::cerr << ss.str();
3707 MInt noTransfer = 0;
3709 if(m_isThermal && !m_isTransport) {
3711 }
else if(m_isTransport && !m_isThermal) {
3713 }
else if(m_isTransport && m_isThermal) {
3721 for(
MInt n = 0; n < noNeighborDomains(); n++) {
3722 sumwin += noWindowCells(n);
3723 sumhalo += noHaloCells(n);
3729 MIntScratchSpace sendNeighOffset(noNeighborDomains(), AT_,
"sendNeighOffset");
3730 MIntScratchSpace recNeighOffset(noNeighborDomains(), AT_,
"recNeighOffset");
3731 sendNeighOffset[0] = 0;
3732 recNeighOffset[0] = 0;
3734 for(
MInt n = 1; n < noNeighborDomains(); n++) {
3735 sendNeighOffset[n] = sendNeighOffset[n - 1] + noWindowCells(n - 1) * noTransfer * m_noDistributions;
3736 recNeighOffset[n] = recNeighOffset[n - 1] + noHaloCells(n - 1) * noTransfer * m_noDistributions;
3739 for(
MInt n = 0; n < noNeighborDomains(); n++) {
3740 MInt sndBuf = sendNeighOffset[n];
3741 for(
MInt j = 0; j < noWindowCells(n); j++) {
3743 MInt offset = j * m_noDistributions +
dist;
3744 sendBuf[sndBuf + offset] = a_oldDistribution(windowCell(n, j),
dist);
3748 sndBuf += m_noDistributions * noWindowCells(n);
3749 for(
MInt j = 0; j < noWindowCells(n); j++) {
3751 MInt offset = j * m_noDistributions +
dist;
3752 sendBuf[sndBuf + offset] = a_oldDistributionThermal(windowCell(n, j),
dist);
3757 sndBuf += m_noDistributions * noWindowCells(n);
3758 for(
MInt j = 0; j < noWindowCells(n); j++) {
3760 MInt offset = j * m_noDistributions +
dist;
3761 sendBuf[sndBuf + offset] = a_oldDistributionTransport(windowCell(n, j),
dist);
3768 for(
MInt n = 0; n < noNeighborDomains(); n++) {
3769 MInt bufsize = noWindowCells(n) * noTransfer * m_noDistributions;
3770 MPI_Issend(&sendBuf[sendNeighOffset[n]], bufsize, MPI_DOUBLE, neighborDomain(n), 0, mpiComm(), &mpi_request[n], AT_,
3771 "sendBuf[sendNeighOffset[n]]");
3777 for(
MInt n = 0; n < noNeighborDomains(); n++) {
3778 MInt bufsize = noHaloCells(n) * noTransfer * m_noDistributions;
3779 MPI_Recv(&recBuf[recNeighOffset[n]], bufsize, MPI_DOUBLE, neighborDomain(n), 0, mpiComm(), &status, AT_,
3780 "recBuf[recNeighOffset[n]]");
3783 for(
MInt n = 0; n < noNeighborDomains(); n++)
3784 MPI_Wait(&mpi_request[n], &status, AT_);
3787 for(
MInt n = 0; n < noNeighborDomains(); n++) {
3788 MInt recBuffer = recNeighOffset[n];
3789 for(
MInt j = 0; j < noHaloCells(n); j++) {
3791 MInt offset = j * m_noDistributions +
dist;
3792 a_oldDistribution(haloCell(n, j),
dist) = recBuf[recBuffer + offset];
3796 recBuffer += m_noDistributions * noHaloCells(n);
3797 for(
MInt j = 0; j < noHaloCells(n); j++) {
3799 MInt offset = j * m_noDistributions +
dist;
3800 a_oldDistributionThermal(haloCell(n, j),
dist) = recBuf[recBuffer + offset];
3805 recBuffer += m_noDistributions * noHaloCells(n);
3806 for(
MInt j = 0; j < noHaloCells(n); j++) {
3808 MInt offset = j * m_noDistributions +
dist;
3809 a_oldDistributionTransport(haloCell(n, j),
dist) = recBuf[recBuffer + offset];
3829 MInt* recalcIdTree) {
3831 if constexpr(nDim != 2 && nDim != 3) {
3832 cerr <<
" In global function saveRestartWithDistributionsPar: wrong number of dimensions !" << endl;
3838 ParallelIo parallelIo(fileName, PIO_REPLACE, mpiComm());
3839 parallelIo.defineScalar(PIO_INT,
"noCells");
3842 const MInt totalNoCells = domainOffset(noDomains());
3845 auto defFloatArray = [&](
const MString arrayName,
const MString varName,
const MInt length) {
3846 parallelIo.defineArray(PIO_FLOAT, arrayName, length);
3847 parallelIo.setAttribute(varName,
"name", arrayName);
3849 auto defIntArray = [&](
const MString arrayName,
const MString varName,
const MInt length) {
3850 parallelIo.defineArray(PIO_INT, arrayName, length);
3851 parallelIo.setAttribute(varName,
"name", arrayName);
3854 auto defineMacroscopicVariables = [&](
const MString name,
const MString prefix,
const MBool old) {
3856 const MString velNames[3] = {
"u",
"v",
"w"};
3857 for(
MInt d = 0; d != nDim; ++d) {
3858 defFloatArray(name + std::to_string(d), prefix + velNames[d], totalNoCells);
3861 defFloatArray(name + std::to_string(nDim), prefix +
"rho", totalNoCells);
3863 if(!m_isTransport && m_isThermal) {
3864 defFloatArray(name + std::to_string(nDim + 1), prefix +
"t", totalNoCells);
3866 if(m_isTransport && !m_isThermal) {
3867 defFloatArray(name + std::to_string(nDim + 1), prefix +
"c", totalNoCells);
3869 if(m_isTransport && m_isThermal) {
3870 defFloatArray(name + std::to_string(nDim + 1), prefix +
"t", totalNoCells);
3871 defFloatArray(name + std::to_string(nDim + 2), prefix +
"c", totalNoCells);
3873 if(m_writeLsData && !old) {
3875 defIntArray(name + std::to_string(nDim + 2), prefix +
"isActive", totalNoCells);
3876 for(
MInt set = 0; set < m_maxNoSets; set++) {
3877 defFloatArray(name + std::to_string(nDim + 3 + set), prefix +
"G_" + std::to_string(set), totalNoCells);
3879 for(
MInt set = 0; set < m_maxNoSets; set++) {
3880 defIntArray(name + std::to_string(nDim + 3 + m_maxNoSets + set), prefix +
"Body_" + std::to_string(set),
3884 defIntArray(name + std::to_string(nDim + 1), prefix +
"isActive", totalNoCells);
3885 for(
MInt set = 0; set < m_maxNoSets; set++) {
3886 defFloatArray(name + std::to_string(nDim + 2 + set), prefix +
"G_" + std::to_string(set), totalNoCells);
3888 for(
MInt set = 0; set < m_maxNoSets; set++) {
3889 defIntArray(name + std::to_string(nDim + 2 + m_maxNoSets + set), prefix +
"Body_" + std::to_string(set),
3895 defFloatArray(name + std::to_string(nDim + 1), prefix +
"alphaG", totalNoCells);
3899 auto defineLocalViscosity = [&](
const MString& name) { defFloatArray(name, name, totalNoCells); };
3902 auto defineDistributions = [&, totalNoCells](
const MString name) {
3903 for(
MInt j = 0; j != m_noDistributions; ++j) {
3904 defFloatArray(name + std::to_string(j), name + std::to_string(j), totalNoCells);
3906 const MString thermalName = name +
"Thermal" + std::to_string(j);
3907 defFloatArray(thermalName, thermalName, totalNoCells);
3910 const MString transportName = name +
"Transport" + std::to_string(j);
3911 defFloatArray(transportName, transportName, totalNoCells);
3916 auto writeMacroscopicVariable = [&](
const MInt index,
const MInt suffix) {
3918 for(
MInt i = 0; i < noInternalCells(); ++i) {
3919 tmp[i] = a_variable(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, index);
3921 parallelIo.writeArray(tmp.getPointer(),
"variables" + std::to_string(suffix));
3924 auto writeMacroscopicOldVariable = [&](
const MInt index,
const MInt suffix) {
3926 for(
MInt i = 0; i < noInternalCells(); ++i) {
3927 tmp[i] = a_oldVariable(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, index);
3929 parallelIo.writeArray(tmp.getPointer(),
"oldVariables" + std::to_string(suffix));
3932 auto writeMacroscopicActiveState = [&](
const MInt suffix) {
3934 for(
MInt i = 0; i < noInternalCells(); ++i) {
3935 tmp[i] = (
MInt)a_isActive(recalcIdTree !=
nullptr ? recalcIdTree[i] : i);
3937 parallelIo.writeArray(tmp.getPointer(),
"variables" + std::to_string(suffix));
3939 auto writeLevelSet = [&](
const MInt set,
const MInt suffix) {
3941 for(
MInt i = 0; i < noInternalCells(); i++) {
3942 tmp[i] = a_levelSetFunctionMB(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, set);
3944 parallelIo.writeArray(tmp.getPointer(),
"variables" + std::to_string(suffix));
3946 auto writeBodyId = [&](
const MInt set,
const MInt suffix) {
3948 for(
MInt i = 0; i < noInternalCells(); i++) {
3949 tmp[i] = a_associatedBodyIds(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, set);
3951 parallelIo.writeArray(tmp.getPointer(),
"variables" + std::to_string(suffix));
3953 auto writeViscosity = [&]() {
3955 for(
MInt i = 0; i < noInternalCells(); ++i) {
3956 tmp[i] = a_nu(recalcIdTree !=
nullptr ? recalcIdTree[i] : i);
3958 parallelIo.writeArray(tmp.getPointer(),
"viscosity");
3962 defineMacroscopicVariables(
"variables",
"",
false);
3963 defineMacroscopicVariables(
"oldVariables",
"old_",
true);
3964 if(m_nonNewtonian) defineLocalViscosity(
"viscosity");
3965 defineDistributions(
"distributions");
3966 defineDistributions(
"oldDistributions");
3967 m_bndCnd->bcDataWriteRestartHeader(parallelIo);
3970 parallelIo.setAttribute(solverId(),
"solverId");
3971 parallelIo.setAttribute(gridInputFileName,
"gridFile",
"");
3972 parallelIo.setAttribute(m_time,
"time");
3974 parallelIo.setAttribute(m_Re,
"ReynoldsNumber");
3975 parallelIo.setAttribute(m_Ma,
"MachNumber");
3977 if(m_controlVelocity) {
3978 parallelIo.setAttribute(m_velocityControl.lastGlobalAvgV,
"lastGlobalAvgV");
3979 parallelIo.setAttribute(m_velocityControl.previousError,
"previousError");
3980 parallelIo.setAttribute(m_velocityControl.integratedError,
"integratedError");
3981 parallelIo.setAttribute(m_velocityControl.derivedError,
"derivedError");
3982 for(
MInt i = 0; i < nDim; i++) {
3983 parallelIo.setAttribute(m_volumeAccel[i],
"volumeAcceleration_" + std::to_string(i));
3988 parallelIo.writeScalar(totalNoCells,
"noCells");
3992 const MPI_Offset firstGlobalId = domainOffset(domainId());
3993 const MPI_Offset localNoCells = noInternalCells();
3994 parallelIo.setOffset(localNoCells, firstGlobalId);
3998 writeMacroscopicVariable(PV->U, suffix++);
3999 writeMacroscopicVariable(PV->V, suffix++);
4000 if constexpr(nDim == 3) writeMacroscopicVariable(PV->W, suffix++);
4001 writeMacroscopicVariable(PV->RHO, suffix++);
4002 if(m_isThermal) writeMacroscopicVariable(PV->T, suffix++);
4003 if(m_isTransport) writeMacroscopicVariable(PV->C, suffix++);
4005 writeMacroscopicActiveState(suffix++);
4006 for(
MInt set = 0; set < m_maxNoSets; set++) {
4007 writeLevelSet(set, suffix++);
4009 for(
MInt set = 0; set < m_maxNoSets; set++) {
4010 writeBodyId(set, suffix++);
4015 for(
MInt i = 0; i < noInternalCells(); ++i) {
4016 const MInt id = recalcIdTree !=
nullptr ? recalcIdTree[i] : i;
4017 tmp_alphaG[i] = a_alphaGas(
id);
4019 parallelIo.writeArray(tmp_alphaG.getPointer(),
"variables" + std::to_string(suffix++));
4021 if(m_nonNewtonian) writeViscosity();
4024 writeMacroscopicOldVariable(PV->U, suffix++);
4025 writeMacroscopicOldVariable(PV->V, suffix++);
4026 if constexpr(nDim == 3) writeMacroscopicOldVariable(PV->W, suffix++);
4027 writeMacroscopicOldVariable(PV->RHO, suffix++);
4028 if(m_isThermal) writeMacroscopicOldVariable(PV->T, suffix++);
4029 if(m_isTransport) writeMacroscopicOldVariable(PV->C, suffix++);
4031 for(
MInt j = 0; j < m_noDistributions; j++) {
4032 const MString distributions =
"distributions" + to_string(j);
4035 for(
MInt i = 0; i < noInternalCells(); ++i)
4036 tmp_distributions[i] = a_distribution(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, j);
4037 parallelIo.writeArray(tmp_distributions.getPointer(), distributions);
4041 for(
MInt j = 0; j < m_noDistributions; j++) {
4042 const MString oldDistributions =
"oldDistributions" + to_string(j);
4044 MFloatScratchSpace tmp_oldDistributions(noInternalCells(), AT_,
"tmp_oldDistributions");
4045 for(
MInt i = 0; i < noInternalCells(); ++i)
4046 tmp_oldDistributions[i] = a_oldDistribution(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, j);
4047 parallelIo.writeArray(tmp_oldDistributions.getPointer(), oldDistributions);
4052 for(
MInt j = 0; j < m_noDistributions; j++) {
4053 const MString distributionsThermal =
"distributionsThermal" + to_string(j);
4055 MFloatScratchSpace tmp_distributionsThermal(noInternalCells(), AT_,
"tmp_distributionsThermal");
4056 for(
MInt i = 0; i < noInternalCells(); ++i)
4057 tmp_distributionsThermal[i] = a_distributionThermal(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, j);
4058 parallelIo.writeArray(tmp_distributionsThermal.getPointer(), distributionsThermal);
4062 for(
MInt j = 0; j < m_noDistributions; j++) {
4063 const MString oldDistributionsThermal =
"oldDistributionsThermal" + to_string(j);
4065 MFloatScratchSpace tmp_oldDistributionsThermal(noInternalCells(), AT_,
"tmp_oldDistributionsThermal");
4066 for(
MInt i = 0; i < noInternalCells(); ++i)
4067 tmp_oldDistributionsThermal[i] = a_oldDistributionThermal(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, j);
4068 parallelIo.writeArray(tmp_oldDistributionsThermal.getPointer(), oldDistributionsThermal);
4074 for(
MInt j = 0; j < m_noDistributions; j++) {
4075 const MString distributionsTransport =
"distributionsTransport" + to_string(j);
4077 MFloatScratchSpace tmp_distributionsTransport(noInternalCells(), AT_,
"tmp_distributionsTransport");
4078 for(
MInt i = 0; i < noInternalCells(); ++i)
4079 tmp_distributionsTransport[i] = a_distributionTransport(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, j);
4080 parallelIo.writeArray(tmp_distributionsTransport.getPointer(), distributionsTransport);
4084 for(
MInt j = 0; j < m_noDistributions; j++) {
4085 const MString oldDistributionsTransport =
"oldDistributionsTransport" + to_string(j);
4087 MFloatScratchSpace tmp_oldDistributionsTransport(noInternalCells(), AT_,
"tmp_oldDistributionsTransport");
4088 for(
MInt i = 0; i < noInternalCells(); ++i)
4089 tmp_oldDistributionsTransport[i] = a_oldDistributionTransport(recalcIdTree !=
nullptr ? recalcIdTree[i] : i, j);
4090 parallelIo.writeArray(tmp_oldDistributionsTransport.getPointer(), oldDistributionsTransport);
4094 m_bndCnd->bcDataWriteRestartData(parallelIo);
4106 stringstream PvFileName;
4107 PvFileName << outputDir() <<
"PV_" << getIdentifier() <<
globalTimeStep << ParallelIo::fileExt();
4109 std::cerr <<
"Writing output file: PV_" << getIdentifier() <<
globalTimeStep <<
" at m_time " << this->m_time
4115 if(m_adaptationSinceLastSolution) {
4116 m_reinitFileName =
"grid" + std::to_string(
globalTimeStep) +
".Netcdf";
4117 m_reinitFilePath = outputDir() + m_reinitFileName;
4118 saveAdaptedGridFile(this->m_recalcIds);
4120 std::vector<MInt> recalcCellIdsSolver(0);
4122 MInt noInternalCellIds;
4123 std::vector<MInt> reorderedCellIds(0);
4124 this->calcRecalcCellIdsSolver(this->m_recalcIds, noCells, noInternalCellIds, recalcCellIdsSolver, reorderedCellIds);
4125 m_adaptationSinceLastSolution =
false;
4126 saveUVWRhoTOnlyPar(PvFileName.str().c_str(), m_reinitFileName.c_str(), recalcCellIdsSolver.data());
4140 for(
MInt i = 0; i < this->maxNoGridCells(); i++) {
4141 p_recalcCellIds[i] = -1;
4145 grid().raw().saveGrid(this->m_reinitFilePath.c_str(), p_recalcCellIds);
4164 if(m_saveDerivatives) {
4169 ASSERT(nDim == 2 || nDim == 3,
"wrong number of dimensions!");
4171 ParallelIo parallelIo(fileName, PIO_REPLACE, mpiComm());
4172 parallelIo.defineScalar(PIO_INT,
"noCells");
4173 parallelIo.setAttribute(solverId(),
"solverId");
4175 const MInt totalNoCells = domainOffset(noDomains());
4178 auto defFloatArray = [&](
const MString& arrayName,
const MString& varName) {
4179 parallelIo.defineArray(PIO_FLOAT, arrayName, totalNoCells);
4180 parallelIo.setAttribute(varName,
"name", arrayName);
4183 auto defineMacroscopicVariables = [&](
const MString& name,
const MString& prefix) {
4186 if constexpr(nDim == 2) {
4187 const MString namesStd[3] = {
"u",
"v",
"rho"};
4189 for(; counter < 3; counter++)
4190 defFloatArray(name + std::to_string(counter), prefix + namesStd[counter]);
4192 const MString namesStd[4] = {
"u",
"v",
"w",
"rho"};
4194 for(; counter < 4; counter++)
4195 defFloatArray(name + std::to_string(counter), prefix + namesStd[counter]);
4199 const MString temperature =
"t";
4200 defFloatArray(name + std::to_string(counter), prefix + temperature);
4205 const MString concentration =
"c";
4206 defFloatArray(name + std::to_string(counter), prefix + concentration);
4209 if(m_particleMomentumCoupling && m_saveExternalForces) {
4210 const MString externalForcing[3] = {
"F_ext_x",
"F_ext_y",
"F_ext_z"};
4211 for(
MInt d = 0; d < nDim; d++, counter++)
4212 defFloatArray(name + std::to_string(counter), prefix + externalForcing[d]);
4216 const MString EELiquid =
"alphaG";
4217 defFloatArray(name + std::to_string(counter), prefix + EELiquid);
4221 if(m_saveDerivatives) {
4222 if constexpr(nDim == 2) {
4224 const MString names[10] = {
"du/dx",
"du/dy",
"dv/dx",
"dv/dy",
"dp_t/dx",
4225 "dp_t/dy",
"vorticity_z",
"Q-criterion",
"Delta-criterion",
"Lambda2"};
4227 for(
MInt d = 0; d < 10; d++, counter++)
4228 defFloatArray(name + std::to_string(counter), prefix + names[d]);
4233 const MString names[19] = {
"du/dx",
"du/dy",
"du/dz",
"dv/dx",
"dv/dy",
4234 "dv/dz",
"dw/dx",
"dw/dy",
"dw/dz",
"dp_t/dx",
4235 "dp_t/dy",
"dp_t/dz",
"vorticity_x",
"vorticity_y",
"vorticity_z",
4236 "vorticity_mag",
"Q-criterion",
"Delta-criterion",
"Lambda2"};
4238 for(
MInt d = 0; d < 19; d++, counter++)
4239 defFloatArray(name + std::to_string(counter), prefix + names[d]);
4242 const MString namesT[3] = {
"dT/dx",
"dT/dy",
"dT/dz"};
4243 for(
MInt d = 0; d < nDim; d++, counter++)
4244 defFloatArray(name + std::to_string(counter), prefix + namesT[d]);
4248 const MString namesC[3] = {
"dC/dx",
"dC/dy",
"dC/dz"};
4249 for(
MInt d = 0; d < nDim; d++, counter++)
4250 defFloatArray(name + std::to_string(counter), prefix + namesC[d]);
4256 defineMacroscopicVariables(
"variables",
"");
4259 MString gridFileName = gridInputFileName;
4260 parallelIo.setAttribute(gridFileName,
"gridFile",
"");
4261 parallelIo.setAttribute(m_time,
"time");
4263 parallelIo.setAttribute(m_Re,
"ReynoldsNumber");
4264 parallelIo.setAttribute(m_Ma,
"MachNumber");
4266 if(m_controlVelocity) {
4267 parallelIo.setAttribute(m_velocityControl.lastGlobalAvgV,
"lastGlobalAvgV");
4268 parallelIo.setAttribute(m_velocityControl.previousError,
"previousError");
4269 parallelIo.setAttribute(m_velocityControl.integratedError,
"integratedError");
4270 parallelIo.setAttribute(m_velocityControl.derivedError,
"derivedError");
4271 for(
MInt i = 0; i < nDim; i++) {
4272 parallelIo.setAttribute(m_volumeAccel[i],
"volumeAcceleration_" + std::to_string(i));
4277 parallelIo.writeScalar(totalNoCells,
"noCells");
4280 const MPI_Offset firstGlobalId = domainOffset(domainId());
4281 const MPI_Offset localNoCells = noInternalCells();
4282 parallelIo.setOffset(localNoCells, firstGlobalId);
4286 MInt start_derivatives = m_noVariables;
4292 for(
MInt d = 0; d != nDim + 1; ++d) {
4294 buffer[
cellId] = a_variable(recalcIdTree !=
nullptr ? recalcIdTree[cellId] : cellId, d);
4296 parallelIo.writeArray(&buffer[0], name + std::to_string(counter++));
4300 if(m_isThermal && !m_isTransport) {
4302 buffer[
cellId] = a_variable(recalcIdTree !=
nullptr ? recalcIdTree[cellId] : cellId, PV->T);
4304 parallelIo.writeArray(&buffer[0], name + std::to_string(counter++));
4306 if(m_isTransport && !m_isThermal) {
4308 buffer[
cellId] = a_variable(recalcIdTree !=
nullptr ? recalcIdTree[cellId] : cellId, PV->C);
4310 parallelIo.writeArray(&buffer[0], name + std::to_string(nDim + 1));
4312 if(m_isTransport && m_isThermal) {
4314 buffer[
cellId] = a_variable(recalcIdTree !=
nullptr ? recalcIdTree[cellId] : cellId, PV->T);
4316 parallelIo.writeArray(&buffer[0], name + std::to_string(nDim + 1));
4319 buffer[
cellId] = a_variable(recalcIdTree !=
nullptr ? recalcIdTree[cellId] : cellId, PV->C);
4321 parallelIo.writeArray(&buffer[0], name + std::to_string(nDim + 2));
4325 if(m_particleMomentumCoupling && m_saveExternalForces) {
4326 for(
MInt d = 0; d < nDim; d++) {
4328 buffer[
cellId] = a_externalForces(recalcIdTree !=
nullptr ? recalcIdTree[cellId] : cellId, d);
4330 parallelIo.writeArray(&buffer[0], name + std::to_string(counter++));
4332 start_derivatives += nDim;
4337 buffer[
cellId] = a_alphaGas(recalcIdTree !=
nullptr ? recalcIdTree[cellId] : cellId);
4339 parallelIo.writeArray(&buffer[0], name + std::to_string(counter++));
4340 start_derivatives += 1;
4345 if(m_saveDerivatives) {
4348 for(
MInt d1 = 0; d1 < nDim; ++d1)
4349 for(
MInt d2 = 0; d2 < nDim; ++d2)
4351 derivatives(d1 * nDim + d2, cellId) = calculateDerivative(cellId, d1, d2);
4353 for(
MInt d = 0; d < nDim; d++) {
4355 if(m_calcTotalPressureGradient != 0) {
4356 derivatives(nDim * nDim + d, cellId) = calculatePressureDerivative(cellId, d);
4358 derivatives(nDim * nDim + d, cellId) = calculateDerivative(cellId, PV->RHO, d);
4364 for(
MInt d1 = 0; d1 < nDim; ++d1)
4365 for(
MInt d2 = 0; d2 < nDim; ++d2)
4366 parallelIo.writeArray(&derivatives(d1 * nDim + d2, 0),
4367 name + std::to_string(start_derivatives + d1 * nDim + d2));
4369 for(
MInt d = 0; d < nDim; d++) {
4370 parallelIo.writeArray(&derivatives(nDim * nDim + d, 0),
4371 name + std::to_string(start_derivatives + nDim * nDim + d));
4380 if constexpr(nDim == 2) {
4381 start_vor = start_derivatives + nDim * nDim + nDim;
4382 start_q = start_vor + 1;
4383 start_d = start_q + 1;
4384 start_lambda2 = start_d + 1;
4386 start_vor = start_derivatives + nDim * nDim + nDim;
4387 start_vor_mag = start_vor + nDim;
4388 start_q = start_vor_mag + 1;
4389 start_d = start_q + 1;
4390 start_lambda2 = start_d + 1;
4395 if constexpr(nDim == 2) {
4398 vorticity_c[
cellId] = derivatives(2, cellId) - derivatives(1, cellId);
4400 parallelIo.writeArray(&vorticity_c[0], name + std::to_string(start_vor));
4402 if constexpr(nDim == 3) {
4405 vorticity_mag[
cellId] = 0.0;
4408 for(
MInt d = 0; d < nDim; ++d) {
4409 MInt v_c1 = (d + 1) % (nDim);
4410 MInt v_c2 = (d + 2) % (nDim);
4413 vorticity_c[
cellId] = derivatives(v_c2 * nDim + v_c1, cellId) - derivatives(v_c1 * nDim + v_c2, cellId);
4417 parallelIo.writeArray(&vorticity_c[0], name + std::to_string(start_vor + d));
4422 vorticity_mag[
cellId] = sqrt(vorticity_mag[cellId]);
4425 parallelIo.writeArray(&vorticity_mag[0], name + std::to_string(start_vor_mag));
4436 getStrainTensor(derivatives, cellId, strain);
4437 getVorticityTensor(derivatives, cellId, vorticity);
4442 q_criterion[
cellId] = 0.5 * (vorticityFrobSq - strainFrobSq);
4445 parallelIo.writeArray(&q_criterion[0], name + std::to_string(start_q));
4452 std::array<std::array<MFloat, nDim>, nDim> dyad;
4454 for(
MInt i = 0; i < nDim; i++) {
4455 for(
MInt j = 0; j < nDim; j++) {
4456 dyad[i][j] = derivatives(j * nDim + i, cellId);
4460 d_criterion[
cellId] = pow((q_criterion[cellId] / 3.0), 3.0) + (0.5 * det) * (0.5 * det);
4462 parallelIo.writeArray(&d_criterion[0], name + std::to_string(start_d));
4473 getStrainTensor(derivatives, cellId, strain);
4474 getVorticityTensor(derivatives, cellId, vorticity);
4481 if constexpr(nDim == 3) {
4483 for(
MInt i = 0; i < nDim; i++) {
4484 for(
MInt j = 0; j < nDim; j++) {
4485 matrix[i][j] = sv(i, j);
4492 for(
MInt i = 0; i < nDim; i++) {
4493 if((evs[i] > evs[((i == 0) ? 0 : (i - 1)) % nDim] && evs[i] < evs[(i + 1) % nDim])
4494 || (evs[i] > evs[(i + 1) % nDim] && evs[i] < evs[((i == 0) ? 0 : (i - 1)) % nDim])) {
4499 lambda_2[
cellId] = evs[l2_index];
4501 TERMM(-1,
"Code is not correct for 2D applications");
4504 parallelIo.writeArray(&lambda_2[0], name + std::to_string(start_lambda2));
4506 if(m_isThermal && !m_isTransport) {
4509 for(
MInt d = 0; d < nDim; d++) {
4511 derivativesT(d, cellId) = calculateDerivative(cellId, PV->T, d);
4513 for(
MInt d = 0; d < nDim; d++) {
4514 parallelIo.writeArray(&derivativesT(d, 0), name + std::to_string(start_lambda2 + 1 + d));
4517 if(m_isTransport && !m_isThermal) {
4520 for(
MInt d = 0; d < nDim; d++) {
4522 derivativesC(d, cellId) = calculateDerivative(cellId, PV->C, d);
4524 for(
MInt d = 0; d < nDim; d++) {
4525 parallelIo.writeArray(&derivativesC(d, 0), name + std::to_string(start_lambda2 + 1 + d));
4528 if(m_isTransport && m_isThermal) {
4531 for(
MInt d = 0; d < nDim; d++) {
4533 derivativesT(d, cellId) = calculateDerivative(cellId, PV->T, d);
4535 for(
MInt d = 0; d < nDim; d++) {
4536 parallelIo.writeArray(&derivativesT(d, 0), name + std::to_string(start_lambda2 + 1 + d));
4541 for(
MInt d = 0; d < nDim; d++) {
4543 derivativesC(d, cellId) = calculateDerivative(cellId, PV->C, d);
4545 for(
MInt d = 0; d < nDim; d++) {
4546 parallelIo.writeArray(&derivativesC(d, 0), name + std::to_string(start_lambda2 + 2 + d));
4563 const MFloat Lb = m_referenceLength;
4566 kpRatio = Context::getSolverProperty<MFloat>(
"kpRatio", this->m_solverId, AT_, &kpRatio);
4568 m_UrmsInit = Context::getSolverProperty<MFloat>(
"UrmsInit", this->m_solverId, AT_);
4569 m_ReLambda = F1 / (F2 * PI * kpRatio) * sqrt(F5 / F6) * m_Re;
4573 m_ReLambda = Context::getSolverProperty<MFloat>(
"ReLambda", this->m_solverId, AT_);
4574 m_UrmsInit = m_ReLambda * 2 * PI * kpRatio * m_nu * sqrt(F6 / F5) / Lb;
4576 mTerm(1, AT_,
"UrmsInit or ReLambda and nu must be specified in property file");
4582 mTerm(1, AT_,
"computeFFTStatistics not implemented for 2D");
4595 constexpr MInt nDim = 3;
4597 std::array<MFloat, nDim * 2> bBox;
4598 m_geometry->getBoundingBox(bBox.data());
4604 const MInt fftLevel = grid().maxUniformRefinementLevel();
4605 if(fftLevel > grid().maxUniformRefinementLevel())
mTerm(1, AT_,
"Non-isotropic mesh is not implemented, yet");
4606 if(fftLevel < minLevel())
mTerm(1, AT_,
"Parents missing for fftLevel < minLevel()");
4608 const MFloat dx = c_cellLengthAtLevel(fftLevel);
4609 const MFloat dxeps = 0.1 * dx;
4610 MInt nx = (bBox[3] - bBox[0] + dxeps) / dx;
4611 MInt ny = (bBox[4] - bBox[1] + dxeps) / dx;
4612 MInt nz = (bBox[5] - bBox[2] + dxeps) / dx;
4623 for(
MInt i = 0; i < noVel; i++) {
4624 pVariables(cellId, i) = a_variable(cellId, i);
4626 oldVariables(cellId, i) = a_oldVariable(cellId, i);
4652 if(a_isHalo(cellId))
continue;
4653 if(a_level(cellId) != fftLevel)
continue;
4663 getReLambdaAndUrmsInit();
4664 if(domainId() == 0) cerr <<
"UrmsInit: " << m_UrmsInit << endl <<
"ReLambdaInit: " << m_ReLambda << endl;
4674 if(a_isHalo(cellId))
continue;
4675 if(a_level(cellId) != fftLevel)
continue;
4677 MFloat actualCellLength = c_cellLengthAtCell(cellId);
4678 MInt xPos = floor((F1B2 * nx + (a_coordinate(cellId, 0) - F1B2 * (actualCellLength)) / dx) + 0.1);
4679 MInt yPos = floor((F1B2 * ny + (a_coordinate(cellId, 1) - F1B2 * (actualCellLength)) / dx) + 0.1);
4680 MInt zPos = floor((F1B2 * nz + (a_coordinate(cellId, 2) - F1B2 * (actualCellLength)) / dx) + 0.1);
4681 if(xPos > nx - 1 || xPos < 0 || yPos > ny - 1 || yPos < 0 || zPos > nz - 1 || zPos < 0) {
4682 cerr <<
"ERROR: wrong array position!" << endl;
4683 cerr <<
"pos=" << xPos <<
", " << yPos <<
", " << zPos << endl;
4684 cerr <<
"coorda = (" << a_coordinate(cellId, 0) <<
", " << a_coordinate(cellId, 1) <<
", "
4685 << a_coordinate(cellId, 2) <<
")" << endl;
4686 cerr <<
"actuallength=" << actualCellLength << endl;
4687 cerr <<
"minlevel=" << minLevel() <<
", maxLevel=" << maxLevel() << endl;
4688 cerr <<
"lenght on level0=" << c_cellLengthAtLevel(0) << endl;
4689 mTerm(1, AT_,
"Wrong array position");
4693 MInt pos = zPos + nz * (yPos + ny * xPos);
4694 for(
MInt i = 0; i < nDim; i++) {
4696 velDat(noLocDat, i) = pVariables(cellId, PV->VV[i]);
4697 velDat(noLocDat, nDim + i) = a_externalForces(cellId, i);
4700 MFloat u1 = pVariables(cellId, PV->VV[i]);
4701 MFloat u0 = oldVariables(cellId, PV->VV[i]);
4702 velDat(noLocDat, 2 * nDim + i) = ((u1 - u0) / F1);
4704 velPos(noLocDat) = pos;
4716 for(
MInt d1 = 0; d1 < nDim; d1++)
4717 strain(d1, d1) = derivatives(d1 * nDim + d1, cellId);
4720 for(
MInt d1 = 0; d1 < nDim - 1; d1++)
4721 for(
MInt d2 = d1 + 1; d2 < nDim; d2++) {
4722 MInt a = d1 * nDim + d2;
4723 MInt b = d2 * nDim + d1;
4724 strain(d1, d2) = 0.5 * (derivatives(
a, cellId) + derivatives(
b, cellId));
4725 strain(d2, d1) = strain(d1, d2);
4736 for(
MInt d1 = 0; d1 < nDim; d1++)
4737 vorticity(d1, d1) = 0.0;
4740 for(
MInt d1 = 0; d1 < nDim - 1; d1++)
4741 for(
MInt d2 = d1 + 1; d2 < nDim; d2++) {
4742 MInt a = d1 * nDim + d2;
4743 MInt b = d2 * nDim + d1;
4744 vorticity(d1, d2) = 0.5 * (derivatives(
a, cellId) - derivatives(
b, cellId));
4745 vorticity(d2, d1) = 0.5 * (derivatives(
b, cellId) - derivatives(
a, cellId));
4769 MInt lbdir1 = 2 * spaceDir;
4770 MInt lbdir2 = lbdir1 + 1;
4772 MInt left = c_neighborId(cellId, lbdir1);
4773 MInt right = c_neighborId(cellId, lbdir2);
4774 MInt leftleft = (c_neighborId(cellId, lbdir1) > -1) ? c_neighborId(left, lbdir1) : -1;
4775 MInt rightright = (c_neighborId(cellId, lbdir2) > -1) ? c_neighborId(right, lbdir2) : -1;
4778 MFloat pt = a_variable(cellId, PV->RHO) * F1B3
4779 + F1B2 / a_variable(cellId, PV->RHO)
4780 * (a_variable(cellId, PV->U) * a_variable(cellId, PV->U)
4781 + a_variable(cellId, PV->V) * a_variable(cellId, PV->V)
4782 + a_variable(cellId, PV->W) * a_variable(cellId, PV->W));
4784 if((left > -1) && (right > -1)) {
4786 a_variable(left, PV->RHO) * F1B3
4787 + F1B2 / a_variable(left, PV->RHO)
4788 * (a_variable(left, PV->U) * a_variable(left, PV->U) + a_variable(left, PV->V) * a_variable(left, PV->V)
4789 + a_variable(left, PV->W) * a_variable(left, PV->W));
4791 MFloat ptRight = a_variable(right, PV->RHO) * F1B3
4792 + F1B2 / a_variable(right, PV->RHO)
4793 * (a_variable(right, PV->U) * a_variable(right, PV->U)
4794 + a_variable(right, PV->V) * a_variable(right, PV->V)
4795 + a_variable(right, PV->W) * a_variable(right, PV->W));
4797 if(leftleft > -1 && rightright > -1) {
4799 MFloat ptLeftLeft = a_variable(leftleft, PV->RHO) * F1B3
4800 + F1B2 / a_variable(leftleft, PV->RHO)
4801 * (a_variable(leftleft, PV->U) * a_variable(leftleft, PV->U)
4802 + a_variable(leftleft, PV->V) * a_variable(leftleft, PV->V)
4803 + a_variable(leftleft, PV->W) * a_variable(leftleft, PV->W));
4805 MFloat ptRightRight = a_variable(rightright, PV->RHO) * F1B3
4806 + F1B2 / a_variable(rightright, PV->RHO)
4807 * (a_variable(rightright, PV->U) * a_variable(rightright, PV->U)
4808 + a_variable(rightright, PV->V) * a_variable(rightright, PV->V)
4809 + a_variable(rightright, PV->W) * a_variable(rightright, PV->W));
4811 gradient = (-ptRightRight + 8.0 * ptRight - 8.0 * ptLeft + ptLeftLeft) / 12.0;
4816 MFloat ptLeftLeft = a_variable(leftleft, PV->RHO) * F1B3
4817 + F1B2 / a_variable(leftleft, PV->RHO)
4818 * (a_variable(leftleft, PV->U) * a_variable(leftleft, PV->U)
4819 + a_variable(leftleft, PV->V) * a_variable(leftleft, PV->V)
4820 + a_variable(leftleft, PV->W) * a_variable(leftleft, PV->W));
4822 gradient = (2.0 * ptRight + 3.0 * pt - 6.0 * ptLeft + ptLeftLeft) / 6.0;
4823 }
else if(rightright > -1) {
4825 MFloat ptRightRight = a_variable(rightright, PV->RHO) * F1B3
4826 + F1B2 / a_variable(rightright, PV->RHO)
4827 * (a_variable(rightright, PV->U) * a_variable(rightright, PV->U)
4828 + a_variable(rightright, PV->V) * a_variable(rightright, PV->V)
4829 + a_variable(rightright, PV->W) * a_variable(rightright, PV->W));
4831 gradient = (-ptRightRight + 6.0 * ptRight - 3.0 * pt - 2.0 * ptLeft) / 6.0;
4833 gradient = (ptRight - ptLeft) / 2.0;
4840 a_variable(left, PV->RHO) * F1B3
4841 + F1B2 / a_variable(left, PV->RHO)
4842 * (a_variable(left, PV->U) * a_variable(left, PV->U) + a_variable(left, PV->V) * a_variable(left, PV->V)
4843 + a_variable(left, PV->W) * a_variable(left, PV->W));
4845 gradient = pt - ptLeft;
4848 else if(right > -1) {
4849 MFloat ptRight = a_variable(right, PV->RHO) * F1B3
4850 + F1B2 / a_variable(right, PV->RHO)
4851 * (a_variable(right, PV->U) * a_variable(right, PV->U)
4852 + a_variable(right, PV->V) * a_variable(right, PV->V)
4853 + a_variable(right, PV->W) * a_variable(right, PV->W));
4855 gradient = ptRight - pt;
4866 TERMM(1,
"Does not work with multi-solver grid concept");
4879 MInt reductionLevel = Context::getSolverProperty<MInt>(
"pp_reductionLevel", m_solverId, AT_);
4883#if defined(MAIA_MS_COMPILER)
4884 MInt recalcIdTreeSize =
4885 grid().raw().reduceToLevel(reductionLevel, [&](
const MInt cellId) { this->interpolateToParentCells(cellId); });
4888 MInt recalcIdTreeSize =
4898 sprintf(buf2,
"%d", reductionLevel);
4900 preName.append(
"Lvl");
4901 preName.append(buf2);
4902 preName.append(
"_");
4904 preName2.append(buf1);
4905 preName2.append(
"_");
4906 preName3 = preName +
"PV_";
4907 preName3.append(buf1);
4908 preName3 = preName3 + ParallelIo::fileExt();
4910 MString tmpFileNameGrid = preName2 + grid().gridInputFileName();
4911 MString tmpFileNameGridPath = outputDir() + tmpFileNameGrid;
4912 MString tmpOutputFileName = outputDir() + preName3;
4916 grid().raw().saveGrid(tmpFileNameGridPath.c_str(), grid().raw().m_haloCells, grid().raw().m_windowCells,
4917 grid().raw().m_azimuthalHaloCells, grid().raw().m_azimuthalUnmappedHaloCells,
4918 grid().raw().m_azimuthalWindowCells, recalcIdTree.begin());
4921 saveUVWRhoTOnlyPar(tmpOutputFileName.c_str(), tmpFileNameGrid.c_str(), recalcIdTree.begin());
4939 vars = &(a_variable(cellId, 0));
4956 if(parentlevel == grid().maxLevel()) {
4957 TERMM(1,
"Interpolation to parent cells residing on the finest level makes no sense, exiting ... ");
4960 for(
MInt id = 0;
id < noInternalCells();
id++) {
4962 if(a_level(
id) != parentlevel || c_isLeafCell(
id))
continue;
4964 for(
MInt v = 0; v < m_noVariables; v++) {
4965 a_variable(
id, v) = .0;
4968 MFloat avg_factor = 1.0 / c_noChildren(
id);
4970 for(
MInt ch = 0; ch <
IPOW2(nDim); ch++) {
4972 if(c_childId(
id, ch) == -1)
continue;
4975 for(
MInt v = 0; v < m_noVariables; v++) {
4976 a_variable(
id, v) += avg_factor * (a_variable(c_childId(
id, ch), v));
4997 if constexpr(nDim != 2 && nDim != 3) {
4998 cerr <<
" In global function loadRestartWithoutDistributionsPar: wrong number of dimensions !" << endl;
5005 ParallelIo parallelIo(fileName, PIO_READ, mpiComm());
5008 parallelIo.getAttribute(&m_time,
"time");
5011 if(m_velocityControl.restart) {
5012 parallelIo.getAttribute(&m_velocityControl.lastGlobalAvgV,
"lastGlobalAvgV");
5013 parallelIo.getAttribute(&m_velocityControl.previousError,
"previousError");
5014 parallelIo.getAttribute(&m_velocityControl.integratedError,
"integratedError");
5015 parallelIo.getAttribute(&m_velocityControl.derivedError,
"derivedError");
5016 for(
MInt i = 0; i < nDim; i++) {
5017 parallelIo.getAttribute(&m_volumeAccel[i],
"volumeAcceleration_" + std::to_string(i));
5022 MPI_Offset dimLen = noInternalCells();
5023 MPI_Offset start = domainOffset(domainId());
5024 parallelIo.setOffset(dimLen, start);
5029 parallelIo.readArray(tmp_velocityU.getPointer(),
"variables0");
5030 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5031 a_variable(i, PV->U) = tmp_velocityU.p[i];
5036 parallelIo.readArray(tmp_velocityV.getPointer(),
"variables1");
5037 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5038 a_variable(i, PV->V) = tmp_velocityV.p[i];
5041 if constexpr(nDim == 3) {
5044 parallelIo.readArray(tmp_velocityW.getPointer(),
"variables2");
5045 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5046 a_variable(i, PV->W) = tmp_velocityW.p[i];
5051 parallelIo.readArray(tmp_rho.getPointer(),
"variables3");
5052 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5053 a_variable(i, PV->RHO) = tmp_rho.p[i];
5056 if(m_isThermal && !m_isTransport) {
5059 parallelIo.readArray(tmp_t.getPointer(),
"variables4");
5060 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5061 a_variable(i, PV->T) = tmp_t.p[i];
5063 if(m_isTransport && !m_isThermal) {
5066 parallelIo.readArray(tmp_c.getPointer(),
"variables4");
5067 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5068 a_variable(i, PV->C) = tmp_c.p[i];
5070 if(m_isTransport && m_isThermal) {
5073 parallelIo.readArray(tmp_t.getPointer(),
"variables4");
5074 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5075 a_variable(i, PV->T) = tmp_t.p[i];
5078 parallelIo.readArray(tmp_c.getPointer(),
"variables5");
5079 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5080 a_variable(i, PV->C) = tmp_c.p[i];
5084 parallelIo.readArray(tmp_alphaG.getPointer(),
"variables4");
5085 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5086 a_alphaGas(i) = tmp_alphaG.p[i];
5091 parallelIo.readArray(tmp_rho.getPointer(),
"variables2");
5092 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5093 a_variable(i, PV->RHO) = tmp_rho.p[i];
5096 if(m_isThermal && !m_isTransport) {
5099 parallelIo.readArray(tmp_t.getPointer(),
"variables3");
5100 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5101 a_variable(i, PV->T) = tmp_t.p[i];
5103 if(m_isTransport && !m_isThermal) {
5106 parallelIo.readArray(tmp_c.getPointer(),
"variables3");
5107 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5108 a_variable(i, PV->C) = tmp_c.p[i];
5110 if(m_isTransport && m_isThermal) {
5113 parallelIo.readArray(tmp_t.getPointer(),
"variables3");
5114 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5115 a_variable(i, PV->T) = tmp_t.p[i];
5118 parallelIo.readArray(tmp_c.getPointer(),
"variables4");
5119 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5120 a_variable(i, PV->C) = tmp_c.p[i];
5126 parallelIo.readArray(tmp_oldVelocityU.getPointer(),
"oldVariables0");
5127 for(
MInt i = 0; i < (
MInt)dimLen; ++i) {
5128 a_oldVariable(i, PV->U) = tmp_oldVelocityU.p[i];
5134 parallelIo.readArray(tmp_oldVelocityV.getPointer(),
"oldVariables1");
5135 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5136 a_oldVariable(i, PV->V) = tmp_oldVelocityV.p[i];
5139 if constexpr(nDim == 3) {
5142 parallelIo.readArray(tmp_oldVelocityW.getPointer(),
"oldVariables2");
5143 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5144 a_oldVariable(i, PV->W) = tmp_oldVelocityW.p[i];
5149 parallelIo.readArray(tmp_oldRho.getPointer(),
"oldVariables3");
5150 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5151 a_oldVariable(i, PV->RHO) = tmp_oldRho.p[i];
5154 if(m_isThermal && !m_isTransport) {
5157 parallelIo.readArray(tmp_oldT.getPointer(),
"oldVariables4");
5158 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5159 a_oldVariable(i, PV->T) = tmp_oldT.p[i];
5162 if(m_isTransport && !m_isThermal) {
5165 parallelIo.readArray(tmp_oldC.getPointer(),
"oldVariables4");
5166 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5167 a_oldVariable(i, PV->C) = tmp_oldC.p[i];
5170 if(m_isTransport && m_isThermal) {
5173 parallelIo.readArray(tmp_oldT.getPointer(),
"oldVariables4");
5174 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5175 a_oldVariable(i, PV->T) = tmp_oldT.p[i];
5179 parallelIo.readArray(tmp_oldC.getPointer(),
"oldVariables5");
5180 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5181 a_oldVariable(i, PV->C) = tmp_oldC.p[i];
5187 parallelIo.readArray(tmp_oldRho.getPointer(),
"oldVariables2");
5188 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5189 a_oldVariable(i, PV->RHO) = tmp_oldRho.p[i];
5192 if(m_isThermal && !m_isTransport) {
5195 parallelIo.readArray(tmp_oldT.getPointer(),
"oldVariables3");
5196 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5197 a_oldVariable(i, PV->T) = tmp_oldT.p[i];
5199 if(m_isTransport && !m_isThermal) {
5202 parallelIo.readArray(tmp_oldC.getPointer(),
"oldVariables3");
5203 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5204 a_oldVariable(i, PV->C) = tmp_oldC.p[i];
5206 if(m_isTransport && m_isThermal) {
5209 parallelIo.readArray(tmp_oldT.getPointer(),
"oldVariables3");
5210 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5211 a_oldVariable(i, PV->T) = tmp_oldT.p[i];
5214 parallelIo.readArray(tmp_oldC.getPointer(),
"oldVariables4");
5215 for(
MInt i = 0; i < (
MInt)dimLen; ++i)
5216 a_oldVariable(i, PV->C) = tmp_oldC.p[i];
5220 m_bndCnd->bcDataReadRestartData(parallelIo);
5243 stringstream fileName;
5246 fileName << outputDir() <<
"restart_" << getIdentifier() <<
globalTimeStep << ParallelIo::fileExt();
5247 grid_ << grid().gridInputFileName();
5249 saveRestartWithDistributionsPar(fileName.str().c_str(), grid().gridInputFileName().c_str());
5260 NEW_TIMER_GROUP(t_restart,
"Restart");
5261 NEW_TIMER(t_restartAll,
"restart file loading", t_restart);
5262 RECORD_TIMER_START(t_restartAll);
5264 m_log <<
"#########################################################################################################"
5267 m_log <<
"## Loading restart file "
5270 m_log <<
"#########################################################################################################"
5276 MBool initializeRestart =
false;
5277 initializeRestart = Context::getSolverProperty<MBool>(
"initRestart", m_solverId, AT_, &initializeRestart);
5278 stringstream varFileName;
5280 m_log <<
" + Restart file information:" << endl;
5281 m_log <<
" - restart time step: " << m_restartTimeStep << endl;
5282 if(m_initFromCoarse)
5283 m_log <<
" - type: "
5284 <<
"from coarse mac. variables " << endl;
5286 m_log <<
" - type: " << (initializeRestart ?
"from mac. variables" :
"from distributions") << endl;
5288 varFileName << restartDir() <<
"restart_" << getIdentifier() << m_restartTimeStep << ParallelIo::fileExt();
5290 m_log <<
" - file name: " << varFileName.str() << endl << endl;
5291 m_log <<
" + Loading the file ..." << endl;
5293 if(m_initFromCoarse)
5294 loadRestartWithoutDistributionsParFromCoarse((varFileName.str()).c_str());
5297 if(!initializeRestart) loadRestartWithDistributionsPar((varFileName.str()).c_str());
5300 loadRestartWithoutDistributionsPar((varFileName.str()).c_str());
5304 if(m_velocityControl.restart) initVolumeForces();
5306 RECORD_TIMER_STOP(t_restartAll);
5328 if(noNeighborDomains() <= 0)
return;
5348 RECORD_TIMER_START(m_t.communication);
5349 (this->*m_receiveMethod)();
5350 RECORD_TIMER_STOP(m_t.communication);
5357 RECORD_TIMER_START(m_t.packing);
5358 (this->*m_gatherMethod)();
5359 RECORD_TIMER_STOP(m_t.packing);
5363 RECORD_TIMER_START(m_t.communication);
5364 (this->*m_sendMethod)();
5365 RECORD_TIMER_STOP(m_t.communication);
5368 RECORD_TIMER_START(m_t.communication);
5369 MPI_Waitall(noNeighborDomains(), &mpi_requestR[0], MPI_STATUS_IGNORE, AT_);
5370 RECORD_TIMER_STOP(m_t.communication);
5373 RECORD_TIMER_START(m_t.unpacking);
5374 (this->*m_scatterMethod)();
5375 RECORD_TIMER_STOP(m_t.unpacking);
5377 RECORD_TIMER_START(m_t.communication);
5378 MPI_Waitall(noNeighborDomains(), &mpi_requestS[0], MPI_STATUS_IGNORE, AT_);
5379 RECORD_TIMER_STOP(m_t.communication);
5400 if(mode == 0)
return;
5403 RECORD_TIMER_START(m_t.packing);
5404 (this->*m_gatherMethod)();
5405 RECORD_TIMER_STOP(m_t.packing);
5407 RECORD_TIMER_START(m_t.communication);
5409 (this->*m_sendMethod)();
5411 (this->*m_receiveMethod)();
5412 RECORD_TIMER_STOP(m_t.communication);
5415 RECORD_TIMER_START(m_t.unpacking);
5416 (this->*m_scatterMethod)();
5417 RECORD_TIMER_STOP(m_t.unpacking);
5423 RECORD_TIMER_START(m_t.solutionStep);
5424 RECORD_TIMER_START(m_t.exchange);
5425 (this->*m_exchangeMethod)(mode);
5426 RECORD_TIMER_STOP(m_t.exchange);
5427 RECORD_TIMER_STOP(m_t.solutionStep);
5449 for(
MInt n = 0; n < noNeighborDomains(); n++) {
5450 for(
MInt var = 0; var < m_noElementsTransfer; var++) {
5451 MInt sendBufferCounter = m_nghbrOffsetsWindow[n][var];
5454 for(
MInt d = 0; d < m_noDistributions; d++) {
5455 MInt id = windowCell(n, j);
5456 MInt offset = sendBufferCounter + (j * m_noDistributions + d);
5457 m_sendBuffers[n][offset] = a_distribution(
id, d);
5460 }
else if(var == 1 && m_isThermal) {
5462 for(
MInt d = 0; d < m_noDistributions; d++) {
5463 MInt id = windowCell(n, j);
5464 MInt offset = sendBufferCounter + (j * m_noDistributions + d);
5465 m_sendBuffers[n][offset] = a_distributionThermal(
id, d);
5468 }
else if(var == 1 && m_isTransport) {
5470 for(
MInt d = 0; d < m_noDistributions; d++) {
5471 MInt id = windowCell(n, j);
5472 MInt offset = sendBufferCounter + (j * m_noDistributions + d);
5473 m_sendBuffers[n][offset] = a_distributionTransport(
id, d);
5476 }
else if(var == 2 && m_isThermal && m_isTransport) {
5478 for(
MInt d = 0; d < m_noDistributions; d++) {
5479 MInt id = windowCell(n, j);
5480 MInt offset = sendBufferCounter + (j * m_noDistributions + d);
5481 m_sendBuffers[n][offset] = a_distributionTransport(
id, d);
5484 }
else if(var == 1 || (var == 2 && (m_isThermal || m_isTransport))
5485 || (var == 3 && m_isThermal && m_isTransport)) {
5487 for(
MInt v = 0; v < m_noVariables; v++) {
5488 MInt id = windowCell(n, j);
5489 MInt offset = sendBufferCounter + (j * m_noVariables + v);
5490 m_sendBuffers[n][offset] = a_variable(
id, v);
5495 for(
MInt v = 0; v < m_noVariables; v++) {
5496 MInt id = windowCell(n, j);
5497 MInt offset = sendBufferCounter + (j * m_noVariables + v);
5498 m_sendBuffers[n][offset] = a_oldVariable(
id, v);
5527 for(
MInt n = 0; n < noNeighborDomains(); n++) {
5528 for(
MInt var = 0; var < m_noElementsTransfer; var++) {
5530 maia::parallelFor<true>(0, noWindowCells(n), [=](
MInt j) {
5531 for(
MInt d = 0; d < (m_noDistributions - 1); d++) {
5532 if(m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] > -1) {
5533 MInt offset = m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d];
5534 m_sendBuffers[n][offset] = a_distribution(windowCell(n, j), d);
5538 }
else if(var == 1 && m_isThermal) {
5539 maia::parallelFor<true>(0, noWindowCells(n), [=](
MInt j) {
5540 for(
MInt d = 0; d < (m_noDistributions - 1); d++) {
5541 if(m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] > -1) {
5543 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] + m_noWindowDistDataPerDomain[n];
5544 m_sendBuffers[n][offset] = a_distributionThermal(windowCell(n, j), d);
5548 }
else if(var == 1 && m_isTransport) {
5549 maia::parallelFor<true>(0, noWindowCells(n), [=](
MInt j) {
5550 for(
MInt d = 0; d < (m_noDistributions - 1); d++) {
5551 if(m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] > -1) {
5553 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] + m_noWindowDistDataPerDomain[n];
5554 m_sendBuffers[n][offset] = a_distributionTransport(windowCell(n, j), d);
5558 }
else if(var == 2 && m_isThermal && m_isTransport) {
5559 maia::parallelFor<true>(0, noWindowCells(n), [=](
MInt j) {
5560 for(
MInt d = 0; d < (m_noDistributions - 1); d++) {
5561 if(m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] > -1) {
5563 m_windowDistsForExchange[n][j * (m_noDistributions - 1) + d] + 2 * m_noWindowDistDataPerDomain[n];
5564 m_sendBuffers[n][offset] = a_distributionTransport(windowCell(n, j), d);
5568 }
else if(var == 1 || (var == 2 && (m_isThermal || m_isTransport))
5569 || (var == 3 && m_isThermal && m_isTransport)) {
5570 maia::parallelFor<true>(0, noWindowCells(n), [=](
MInt j) {
5571 for(
MInt v = 0; v < m_noVariables; v++) {
5572 MInt id = windowCell(n, j);
5573 MInt offset = (j * m_noVariables + v) + (m_noDistsTransfer * m_noWindowDistDataPerDomain[n]);
5574 m_sendBuffers[n][offset] = a_variable(
id, v);
5578 maia::parallelFor<true>(0, noWindowCells(n), [=](
MInt j) {
5579 for(
MInt v = 0; v < m_noVariables; v++) {
5580 MInt id = windowCell(n, j);
5581 MInt offset = (j * m_noVariables + v) + (m_noDistsTransfer * m_noWindowDistDataPerDomain[n])
5582 + (noWindowCells(n) * m_noVariables);
5583 m_sendBuffers[n][offset] = a_oldVariable(
id, v);
5611 for(
MInt n = 0; n < noNeighborDomains(); n++) {
5612 for(
MInt var = 0; var < m_noElementsTransfer; var++) {
5613 MInt receiveBufferCounter = m_nghbrOffsetsHalo[n][var];
5616 for(
MInt d = 0; d < m_noDistributions; d++) {
5617 MInt id = haloCell(n, j);
5618 MInt offset = receiveBufferCounter + (j * m_noDistributions + d);
5619 a_distribution(
id, d) = m_receiveBuffers[n][offset];
5622 }
else if(var == 1 && m_isThermal) {
5624 for(
MInt d = 0; d < m_noDistributions; d++) {
5625 MInt id = haloCell(n, j);
5626 MInt offset = receiveBufferCounter + (j * m_noDistributions + d);
5627 a_distributionThermal(
id, d) = m_receiveBuffers[n][offset];
5630 }
else if(var == 1 && m_isTransport) {
5632 for(
MInt d = 0; d < m_noDistributions; d++) {
5633 MInt id = haloCell(n, j);
5634 MInt offset = receiveBufferCounter + (j * m_noDistributions + d);
5635 a_distributionTransport(
id, d) = m_receiveBuffers[n][offset];
5638 }
else if(var == 2 && m_isThermal && m_isTransport) {
5640 for(
MInt d = 0; d < m_noDistributions; d++) {
5641 MInt id = haloCell(n, j);
5642 MInt offset = receiveBufferCounter + (j * m_noDistributions + d);
5643 a_distributionTransport(
id, d) = m_receiveBuffers[n][offset];
5646 }
else if(var == 1 || (var == 2 && (m_isThermal || m_isTransport))
5647 || (var == 3 && m_isThermal && m_isTransport)) {
5649 for(
MInt v = 0; v < m_noVariables; v++) {
5650 MInt id = haloCell(n, j);
5651 MInt offset = receiveBufferCounter + (j * m_noVariables + v);
5652 a_variable(
id, v) = m_receiveBuffers[n][offset];
5657 for(
MInt v = 0; v < m_noVariables; v++) {
5658 MInt id = haloCell(n, j);
5659 MInt offset = receiveBufferCounter + (j * m_noVariables + v);
5660 a_oldVariable(
id, v) = m_receiveBuffers[n][offset];
5689 for(
MInt n = 0; n < noNeighborDomains(); n++) {
5690 for(
MInt var = 0; var < m_noElementsTransfer; var++) {
5692 maia::parallelFor<true>(0, noHaloCells(n), [=](
MInt j) {
5693 for(
MInt d = 0; d < (m_noDistributions - 1); d++) {
5694 if(m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] > -1) {
5695 MInt offset = m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d];
5696 a_distribution(haloCell(n, j), d) = m_receiveBuffers[n][offset];
5700 }
else if(var == 1 && m_isThermal) {
5701 maia::parallelFor<true>(0, noHaloCells(n), [=](
MInt j) {
5702 for(
MInt d = 0; d < (m_noDistributions - 1); d++) {
5703 if(m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] > -1) {
5704 MInt offset = m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] + m_noHaloDistDataPerDomain[n];
5705 a_distributionThermal(haloCell(n, j), d) = m_receiveBuffers[n][offset];
5709 }
else if(var == 1 && m_isTransport) {
5710 maia::parallelFor<true>(0, noHaloCells(n), [=](
MInt j) {
5711 for(
MInt d = 0; d < (m_noDistributions - 1); d++) {
5712 if(m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] > -1) {
5713 MInt offset = m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] + m_noHaloDistDataPerDomain[n];
5714 a_distributionTransport(haloCell(n, j), d) = m_receiveBuffers[n][offset];
5718 }
else if(var == 2 && m_isThermal && m_isTransport) {
5719 maia::parallelFor<true>(0, noHaloCells(n), [=](
MInt j) {
5720 for(
MInt d = 0; d < (m_noDistributions - 1); d++) {
5721 if(m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] > -1) {
5723 m_haloDistsForExchange[n][j * (m_noDistributions - 1) + d] + 2 * m_noHaloDistDataPerDomain[n];
5724 a_distributionTransport(haloCell(n, j), d) = m_receiveBuffers[n][offset];
5728 }
else if(var == 1 || (var == 2 && (m_isThermal || m_isTransport))
5729 || (var == 3 && m_isThermal && m_isTransport)) {
5730 maia::parallelFor<true>(0, noHaloCells(n), [=](
MInt j) {
5731 for(
MInt v = 0; v < m_noVariables; v++) {
5732 MInt id = haloCell(n, j);
5733 MInt offset = (j * m_noVariables + v) + m_noDistsTransfer * m_noHaloDistDataPerDomain[n];
5734 a_variable(
id, v) = m_receiveBuffers[n][offset];
5738 maia::parallelFor<true>(0, noHaloCells(n), [=](
MInt j) {
5739 for(
MInt v = 0; v < m_noVariables; v++) {
5740 MInt id = haloCell(n, j);
5741 MInt offset = (j * m_noVariables + v) + m_noDistsTransfer * m_noHaloDistDataPerDomain[n]
5742 + (noHaloCells(n) * m_noVariables);
5743 a_oldVariable(
id, v) = m_receiveBuffers[n][offset];
5767 for(
MInt i = 0; i < noNeighborDomains(); i++) {
5768 bufSize = noWindowCells(i) * m_dataBlockSizeTotal;
5769 MPI_Issend(m_sendBuffers[i], bufSize, MPI_DOUBLE, neighborDomain(i), 0, mpiComm(), &mpi_request[i], AT_,
5770 "m_sendBuffers[i]");
5786 for(
MInt i = 0; i < noNeighborDomains(); i++) {
5788 m_noDistsTransfer * m_noWindowDistDataPerDomain[i] + noWindowCells(i) * m_noVariables * m_noVarsTransfer;
5789 if(totalData == 0)
continue;
5790 MPI_Issend(m_sendBuffers[i], totalData, MPI_DOUBLE, neighborDomain(i), 0, mpiComm(), &mpi_request[i], AT_,
5791 "m_sendBuffers[i]");
5810 for(
MInt i = 0; i < noNeighborDomains(); i++) {
5811 const MInt bufSize = noHaloCells(i) * m_dataBlockSizeTotal;
5812 MPI_Recv(m_receiveBuffers[i], bufSize, MPI_DOUBLE, neighborDomain(i), 0, mpiComm(), &status, AT_,
5813 "m_receiveBuffers[i]");
5815 for(
MInt i = 0; i < noNeighborDomains(); i++)
5816 MPI_Wait(&mpi_request[i], &status, AT_);
5834 for(
MInt i = 0; i < noNeighborDomains(); i++) {
5836 m_noDistsTransfer * m_noHaloDistDataPerDomain[i] + noHaloCells(i) * m_noVariables * m_noVarsTransfer;
5837 if(totalData == 0)
continue;
5838 MPI_Recv(m_receiveBuffers[i], totalData, MPI_DOUBLE, neighborDomain(i), 0, mpiComm(), &status, AT_,
5839 "m_receiveBuffers[i]");
5841 for(
MInt i = 0; i < noNeighborDomains(); i++) {
5843 m_noDistsTransfer * m_noHaloDistDataPerDomain[i] + noHaloCells(i) * m_noVariables * m_noVarsTransfer;
5844 if(totalData == 0)
continue;
5845 MPI_Wait(&mpi_request[i], MPI_STATUS_IGNORE, AT_);
5867 m_log <<
" + Setting active cell list" << endl;
5869 setIsActiveDefaultStates();
5870 fillActiveCellList();
5882 m_activeCellList =
nullptr;
5883 m_currentMaxNoCells = (m_cells.size()
5898 m_activeCellList =
nullptr;
5900 m_currentMaxNoCells = 0;
5916 m_currentMaxNoCells = m_cells.size();
5924 NEW_SUB_TIMER(t_ifaceCells,
"interface cells", m_t.solver);
5925 RECORD_TIMER_START(t_ifaceCells);
5928 m_log <<
" + Interface cell generation" << endl;
5930 resetInterfaceCells();
5931 initializeInterfaceCells();
5932 buildInterfaceCells();
5933 setInterpolationNeighbors();
5934 setInterpolationCoefficients();
5936 RECORD_TIMER_STOP(t_ifaceCells);
5942 if(m_isRefined && m_correctInterfaceBcCells) {
5943 setInterpolationNeighborsBC();
5944 setInterpolationCoefficientsBC();
5957 const MInt interfaceLevels = maxLevel() - minLevel();
5959 if(interfaceLevels == 0)
return;
5971 m_interfaceCellSize = 0.5;
5972 m_interfaceCellSize = Context::getSolverProperty<MFloat>(
"interfaceCellSize", m_solverId, AT_, &m_interfaceCellSize);
5974 if(m_interfaceCellSize > noInternalCells())
5976 " LbSolver::initializeInterfaceCells: You have requested more space for interface cells than for "
5978 "internal cells. What is the sense of that? Exiting!");
5980 MInt nointerfaceCellSize = (
MInt)(m_interfaceCellSize * noInternalCells());
5982 m_log <<
" - initializing interface cells for levels " << minLevel() <<
" to " << maxLevel() << endl;
5983 m_log <<
" * currently using a size of " << nointerfaceCellSize <<
" for interface child and parent cells"
5985 m_log <<
" * ratio: " << m_interfaceCellSize << endl;
5986 m_log <<
" * collector size: "
5993 for(
MInt l = 0; l < interfaceLevels; l++) {
5994 m_interfaceChildren.emplace_back();
5997 for(
MInt i = 0; i < interfaceLevels; i++) {
6002 for(
MInt l = 0; l < interfaceLevels; l++) {
6003 m_interfaceParents.emplace_back();
6006 for(
MInt i = 0; i < interfaceLevels; i++) {
6018 m_log <<
" - resetting interface cells" << endl;
6019 const MInt interfaceLevels = maxLevel() - minLevel();
6021 if(interfaceLevels == 0)
return;
6023 MInt noCells = m_cells.size();
6027 for(
MInt i = 0; i < noCells; i++) {
6028 a_isInterfaceChild(i) =
false;
6029 a_isInterfaceParent(i) =
false;
6035 if(m_interfaceChildren.size() != 0) {
6036 for(
auto&& children : m_interfaceChildren) {
6039 m_interfaceChildren.clear();
6041 if(m_interfaceParents.size() != 0) {
6042 for(
auto&& parents : m_interfaceParents) {
6045 m_interfaceParents.clear();
6058 m_log <<
" - building interface cells" << endl;
6060 const MInt interfaceLevels = maxLevel() - minLevel();
6062 if(interfaceLevels == 0)
return;
6064 MInt noCells = m_cells.size();
6066 MInt pointedId = -1;
6072 if constexpr(nDim == 3) {
6080 m_log <<
" * finding interface children" << endl;
6082 mAlloc(m_noInterfaceChildren, interfaceLevels,
"m_noInterfaceChildren", 0, AT_);
6085 for(
MInt i = 0; i < noInternalCells(); i++) {
6087 if(a_level(i) == minLevel())
continue;
6094 if(c_parentId(i) > -1) {
6095 parentId = c_parentId(i);
6098 const MInt currentInterfaceLevel = a_level(pointedId) - minLevel() - 1;
6102 if(a_hasNeighbor(pointedId,
dist) == 0 && a_hasNeighbor(parentId,
dist)) {
6104 if(c_isLeafCell(c_neighborId(parentId,
dist))) {
6107 a_isInterfaceChild(pointedId) =
true;
6108 m_noInterfaceChildren[currentInterfaceLevel]++;
6114 if(!m_correctInterfaceBcCells && a_isInterfaceChild(pointedId)) {
6117 if(a_hasNeighbor(pointedId,
dist) == 0 && a_hasNeighbor(parentId,
dist) == 0) {
6118 a_isInterfaceChild(pointedId) =
false;
6119 m_noInterfaceChildren[currentInterfaceLevel]--;
6127 for(
MInt level = 0; level < interfaceLevels; level++)
6128 if(m_noInterfaceChildren[level] > 0)
6129 m_log <<
" # found children for level " << minLevel() + level + 1 <<
": " << m_noInterfaceChildren[level]
6135 m_log <<
" # writing children to collector " << endl;
6137 for(
MInt i = 0; i < noCells; i++) {
6140 if(!a_isInterfaceChild(i))
continue;
6142 parentId = c_parentId(i);
6146 const MInt currentInterfaceLevel = a_level(i) - minLevel() - 1;
6147 m_interfaceChildren[currentInterfaceLevel]->append();
6149 auto& currentInterfaceCell =
6150 m_interfaceChildren[currentInterfaceLevel]->a[m_interfaceChildren[currentInterfaceLevel]->size() - 1];
6151 currentInterfaceCell.m_cellId = i;
6154 currentInterfaceCell.m_position = 0;
6155 for(
MInt dim = 0; dim < nDim; dim++) {
6156 if(a_coordinate(pointedId, dim) < a_coordinate(parentId, dim)) {
6157 currentInterfaceCell.m_position += 0;
6159 currentInterfaceCell.m_position +=
IPOW2(dim);
6166 m_log <<
" * finding interface parents" << endl;
6168 mAlloc(m_noInterfaceParents, interfaceLevels,
"m_noInterfaceParents", 0, AT_);
6171 for(
MInt i = 0; i < noCells; i++) {
6174 if(a_level(i) == maxLevel())
continue;
6176 if(c_isLeafCell(i))
continue;
6179 if(c_childId(i, k) < 0 || c_childId(i, k) >= noCells) {
6180 if(c_childId(i, k) != -1) cout <<
"STRANGE ID " << c_childId(i, k) << endl;
6183 if(a_isInterfaceChild(c_childId(i, k))) {
6186 if(c_noChildren(i) !=
IPOW2(nDim) && i < noInternalCells()) {
6187 stringstream errorMessage;
6188 errorMessage <<
" LbInterface::identifyInterfaceCells: Insufficient number of children in interface "
6189 "parent cell with id "
6190 << i <<
", Exiting!";
6191 TERMM(1, errorMessage.str());
6194 a_isInterfaceParent(i) =
true;
6196 const MInt currentInterfaceLevel = a_level(i) - minLevel();
6197 m_noInterfaceParents[currentInterfaceLevel]++;
6204 for(
MInt level = 0; level < interfaceLevels; level++)
6205 if(m_noInterfaceParents[level] > 0)
6206 m_log <<
" # found parents for level " << minLevel() + level + 1 <<
": " << m_noInterfaceParents[level]
6210 m_log <<
" * adding halo cells " << endl;
6213 for(
MInt i = noInternalCells(); i < noCells; i++) {
6215 if(a_level(i) == maxLevel())
continue;
6218 if(c_isLeafCell(i))
continue;
6221 emptyParent =
false;
6222 for(
MInt child = 0; child <
IPOW2(nDim); child++) {
6223 if(c_childId(i, child) < 0) {
6228 if(emptyParent)
continue;
6230 for(
MInt k = 0; k < maxNoNghbrs; k++) {
6231 if(a_hasNeighbor(i, k) > 0 && c_isLeafCell(c_neighborId(i, k))) {
6232 a_isInterfaceParent(i) =
true;
6234 const MInt currentInterfaceLevel = a_level(i) - minLevel();
6235 m_noInterfaceParents[currentInterfaceLevel]++;
6241 for(
MInt level = 0; level < interfaceLevels; level++)
6242 if(m_noInterfaceParents[level] > 0)
6243 m_log <<
" # new number of parents for level " << minLevel() + level + 1 <<
": "
6244 << m_noInterfaceParents[level] << endl;
6246 m_log <<
" # writing parents to collector " << endl;
6250 for(
MInt i = 0; i < noCells; i++) {
6254 if(!a_isInterfaceParent(i))
continue;
6257 const MInt currentInterfaceLevel = a_level(i) - minLevel();
6258 m_interfaceParents[currentInterfaceLevel]->append();
6260 auto& currentParentCell =
6261 m_interfaceParents[currentInterfaceLevel]->a[m_interfaceParents[currentInterfaceLevel]->size() - 1];
6262 currentParentCell.m_cellId = i;
6274 m_log <<
" - setting interpolation coefficients" << endl;
6276 MInt noInterpolationNghbrs;
6277 noInterpolationNghbrs =
IPOW2(nDim);
6280 switch(m_interpolationType) {
6283 for(
MInt j = 0; j < noInterpolationNghbrs; j++) {
6292 " In function LbInterfaceD2Q9::setInterpolationCoefficients() : quadratic interpolation is not "
6293 "implemented yet!");
6299 " In function LbInterfaceD2Q9::setInterpolationCoefficients() : cubic interpolation is not "
6300 "implemented yet!");
6304 TERMM(1,
" In function LbInterfaceD2Q9::setInterpolationCoefficients() : Unknown interpolation type!");
6309 for(
MInt i = 0; i < (maxLevel() - minLevel()); i++) {
6310 for(
MInt j = 0; j < m_interfaceChildren[i]->size(); j++) {
6311 position = m_interfaceChildren[i]->a[j].m_position;
6312 for(
MInt k = 0; k < noInterpolationNghbrs; k++) {
6313 m_interfaceChildren[i]->a[j].m_interpolationCoefficients[k] = interpolationCoefficients(position, k);
6327 m_log <<
" - setting interpolation neighbors" << endl;
6333 if constexpr(nDim == 2) {
6335 switch(m_interpolationType) {
6337 for(
MInt i = 0; i < maxLevel() - minLevel(); i++) {
6338 for(
MInt j = 0; j < m_noInterfaceChildren[i]; j++) {
6339 currentId = m_interfaceChildren[i]->a[j].m_cellId;
6340 position = m_interfaceChildren[i]->a[j].m_position;
6343 if(dir < IPOW3[nDim] - 1) {
6346 if(a_hasNeighbor(currentId, dir) == 0) {
6347 m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k] = c_neighborId(c_parentId(currentId), dir);
6349 m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k] = c_parentId(c_neighborId(currentId, dir));
6352 m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k] = c_parentId(currentId);
6362 " In function LbInterfaceD2Q9::setInterpolationCoefficients() : quadratic interpolation not "
6363 "implemented yet!");
6368 " In function LbInterfaceD2Q9::setInterpolationCoefficients() : cubic interpolation not "
6369 "implemented yet!");
6372 TERMM(1,
" In function LbInterfaceD2Q9::setInterpolationCoefficients() : Unknown interpolation type!");
6379 for(
MInt i = 0; i < maxLevel() - minLevel(); i++) {
6380 for(
MInt j = 0; j < m_noInterfaceChildren[i]; j++) {
6381 currentId = m_interfaceChildren[i]->a[j].m_cellId;
6382 position = m_interfaceChildren[i]->a[j].m_position;
6385 if(dir < IPOW3[nDim] - 1) {
6387 const MInt parentId = c_parentId(currentId);
6388 const MInt nghbrId = c_neighborId(parentId, dir);
6389 m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k] = nghbrId;
6392 m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k] = c_parentId(currentId);
6409 std::cout <<
"WARNING: Refined boundary cells are problematic in 2D! "
6410 <<
"This is not implemented, yet." << std::endl;
6414 for(
MInt i = 0; i < (maxLevel() - minLevel()); i++) {
6416 for(
MInt j = 0; j < m_noInterfaceChildren[i]; j++) {
6417 const MInt cellId = m_interfaceChildren[i]->a[j].m_cellId;
6424 MBool performCorrection =
false;
6425 constexpr MInt maxNoNghbrs = IPOW3[nDim] - 1;
6426 std::array<MBool, maxNoNghbrs> validShiftDirection;
6427 validShiftDirection.fill(
true);
6428 for(
MInt k = 0; k < m_interfaceChildren[i]->a[j].m_noInterpolationNeighbors; k++) {
6429 const MInt nghbrId = m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k];
6430 if(nghbrId < 0 || !a_isActive(nghbrId) || a_isBndryCell(nghbrId)) {
6431 performCorrection =
true;
6434 const MInt srcId = (nghbrId < 0) ? c_parentId(cellId) : nghbrId;
6436 for(
MInt d = 0; d < maxNoNghbrs; d++) {
6437 if(a_hasNeighbor(srcId, d)) {
6438 const MInt srcNghbrId = c_neighborId(srcId, d);
6439 if(!a_isActive(srcNghbrId) || a_isBndryCell(srcNghbrId)) {
6440 validShiftDirection[d] =
false;
6443 validShiftDirection[d] =
false;
6449 if(performCorrection) {
6450 MBool validShift =
false;
6451 for(
MInt d = 0; d < maxNoNghbrs; d++) {
6452 if(validShiftDirection[d]) {
6454 for(
MInt k = 0; k < m_interfaceChildren[i]->a[j].m_noInterpolationNeighbors; k++) {
6455 const MInt nghbrId = m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k];
6457 m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k] = c_neighborId(nghbrId, d);
6459 std::array<MInt, nDim> newDirVec;
6460 const MInt position = m_interfaceChildren[i]->a[j].m_position;
6461 for(
MInt l = 0; l < nDim; l++) {
6470 const MInt parentId = c_parentId(cellId);
6471 if(newDir == IPOW3[nDim] - 1) {
6472 m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k] = parentId;
6474 m_interfaceChildren[i]->a[j].m_interpolationNeighbors[k] = c_neighborId(parentId, newDir);
6482 std::stringstream err;
6483 err <<
"ERROR: interpolation stencil for BC interface could not be corrected!";
6484 err <<
" [globalId: " << c_globalId(cellId) <<
"]" << std::endl;
6485 TERMM(1, err.str());
6488 InterfaceLink tmp = {i, j};
6489 m_interfaceChildrenBc.push_back(tmp);
6497 TERMM(1,
"Only nDim=2 and nDim=3 allowed!");
6517 for(
auto link : m_interfaceChildrenBc) {
6518 const MInt i = link.levelId;
6519 const MInt j = link.interfaceId;
6523 std::array<MFloat, nDim> c0;
6524 const MInt cellId = m_interfaceChildren[i]->a[j].m_cellId;
6525 const MInt cellId0 = m_interfaceChildren[i]->a[j].m_interpolationNeighbors[0];
6526 for(
MInt d = 0; d < nDim; d++) {
6527 const MInt cellId1 = m_interfaceChildren[i]->a[j].m_interpolationNeighbors[
IPOW2(d)];
6528 const MFloat x = a_coordinate(cellId, d);
6529 const MFloat x0 = a_coordinate(cellId0, d);
6530 const MFloat x1 = a_coordinate(cellId1, d);
6531 c0[d] = (x1 - x) / (x1 - x0);
6534 for(
MInt m = 0; m < 8; m++) {
6537 for(
MInt d = nDim - 1; d > -1; d--) {
6539 tmp = tmp %
IPOW2(d);
6540 const MBool chk = (chk0 == 0);
6541 weight *= (chk) ? c0[d] : 1 - c0[d];
6543 m_interfaceChildren[i]->a[j].m_interpolationCoefficients[m] = weight;
6549 TERMM(1,
"Only nDim=2 and nDim=3 allowed!");
6566 MString filename = outputDir() +
"PV_";
6568 sprintf(buf,
"%d", timeStep);
6569 filename.append(buf);
6570 filename += ParallelIo::fileExt();
6572 loadRestartWithoutDistributionsPar(filename.c_str());
6588#ifdef WAR_NVHPC_PSTL
6590 TERMM(1,
"getSampleVariables: var runs out of scope after this function call!");
6592 MFloat var[nDim + 2] = {F0};
6593 for(
MInt d = 0; d < nDim + 2; d++) {
6594 var[d] = a_variable(cellId, d);
6598 MFloat var[nDim + 1] = {F0};
6599 for(
MInt d = 0; d < nDim + 1; d++) {
6600 var[d] = a_variable(cellId, d);
6605 vars = a_variables_ptr(cellId);
6611 const MInt noVars = vars.size();
6612 ASSERT(noVars <= m_noVariables, "noVars > m_noVariables
");
6613 for(MInt varId = 0; varId < noVars; varId++) {
6614 vars[varId] = a_variable(cellId, varId);
6619MBool LbSolver<nDim>::getSampleVarsDerivatives(const MInt cellId, std::vector<MFloat>& vars) {
6620 auto isValid = [&](const MInt l_cellId, const MInt l_dir) -> MInt {
6621 if(l_cellId != -1) {
6622 const MInt nghbrId = c_neighborId(l_cellId, l_dir);
6623 if(nghbrId > -1 && a_isActive(nghbrId)) {
6629 const MFloat F1bdx = FFPOW2(maxLevel() - a_level(cellId)); // LB units
6630 for(MInt dir = 0; dir < nDim; dir++) {
6631 const MInt dir0 = 2 * dir;
6632 const MInt dir1 = 2 * dir + 1;
6633 // 1) set coefficients
6634 constexpr MInt noCellsInStencil = 5;
6635 std::array<MFloat, noCellsInStencil> coef;
6636 std::array<MInt, noCellsInStencil> cellIds;
6637 // 1.1) get cellIds and coefficients
6638 cellIds[1] = isValid(cellId, dir0); // left
6639 cellIds[0] = isValid(cellIds[1], dir0); // leftleft
6640 cellIds[2] = cellId; // middle
6641 cellIds[3] = isValid(cellId, dir1); // right
6642 cellIds[4] = isValid(cellIds[3], dir1); // rightright
6643 getDerivativeStencilAndCoefficient(cellIds, coef);
6644 // 1.2) scale coefficients with correct dx in case of refined
6645 for(MInt i = 0; i < noCellsInStencil; i++) {
6648 // 2) calculate gradients
6650 for(MInt veloId = 0; veloId < nDim; veloId++) {
6651 const MInt index = nDim * veloId + dir;
6653 for(MInt i = 0; i < noCellsInStencil; i++) {
6654 vars[index] += cellIds[i] < 0 ? 0.0 : coef[i] * a_variable(cellIds[i], veloId);
6658 for(MInt varId = nDim; varId < m_noVariables; varId++) {
6659 const MInt index = nDim * varId + dir;
6661 for(MInt i = 0; i < noCellsInStencil; i++) {
6662 vars[index] += cellIds[i] < 0 ? 0.0 : coef[i] * a_variable(cellIds[i], varId);
6667 const MInt index = m_noVariables * nDim + dir;
6668 for(MInt i = 0; i < noCellsInStencil; i++) {
6669 vars[index] += cellIds[i] < 0 ? 0.0 : coef[i] * a_alphaGas(cellIds[i]);
6677void LbSolver<nDim>::storeOldDistributions() {
6678 if(!m_cells.savePrevVars()) {
6679 TERMM(1, "Cell collector not configured to stor previous values!
");
6681 for(MInt cellId = 0; cellId < a_noCells(); cellId++) {
6682 for(MInt distr = 0; distr < m_noDistributions; distr++) {
6683 a_previousDistribution(cellId, distr) = a_oldDistribution(cellId, distr);
6689void LbSolver<nDim>::storeOldVariables() {
6690 if(!m_cells.savePrevVars()) {
6691 TERMM(1, "Cell collector not configured to stor previous values!
");
6693 for(MInt cellId = 0; cellId < a_noCells(); cellId++) {
6694 a_previousVariable(cellId, PV->RHO) = a_oldVariable(cellId, PV->RHO);
6695 a_previousVariable(cellId, PV->U) = a_oldVariable(cellId, PV->U);
6696 a_previousVariable(cellId, PV->V) = a_oldVariable(cellId, PV->V);
6697 a_previousVariable(cellId, PV->W) = a_oldVariable(cellId, PV->W);
6702void LbSolver<nDim>::initNu(const MInt cellId, const MFloat nu) {
6704 if(m_cells.saveNuT()) {
6705 a_nuT(cellId) = 0.0;
6706 if(m_cells.saveOldNu()) {
6707 a_oldNuT(cellId) = 0.0;
6710 if(m_cells.saveOldNu()) {
6711 a_oldNu(cellId) = nu;
6717MBool LbSolver<nDim>::solutionStep() {
6720 RECORD_TIMER_START(m_t.solutionStep);
6722 // calls the exchange method
6723 RECORD_TIMER_START(m_t.exchange);
6724 (this->*m_exchangeMethod)(0);
6725 RECORD_TIMER_STOP(m_t.exchange);
6727 // update macroscopic variables
6728 if(m_updateMacroscopicLocation == PRECOLLISION) {
6729 updateVariablesFromOldDist_preCollision();
6732 // update/reset viscosity
6735 // pre-collison source term call
6736 RECORD_TIMER_START(m_t.srcTerms);
6737 preCollisionSrcTerm();
6738 RECORD_TIMER_STOP(m_t.srcTerms);
6741 RECORD_TIMER_START(m_t.collision);
6742 (this->*m_solutionStepMethod)();
6743 RECORD_TIMER_STOP(m_t.collision);
6745 // post-collison source term call
6746 RECORD_TIMER_START(m_t.srcTerms);
6747 postCollisionSrcTerm();
6748 RECORD_TIMER_STOP(m_t.srcTerms);
6750 // post-collision boundary conditions
6751 RECORD_TIMER_START(m_t.collisionBC);
6753 RECORD_TIMER_STOP(m_t.collisionBC);
6755 // calls the exchange method
6756 RECORD_TIMER_START(m_t.exchange);
6757 (this->*m_exchangeMethod)(1);
6758 RECORD_TIMER_STOP(m_t.exchange);
6761 RECORD_TIMER_START(m_t.propagation);
6762 (this->*m_propagationStepMethod)();
6763 RECORD_TIMER_STOP(m_t.propagation);
6765 // TODO labels:LB @johannes/julian: I made this step in updateVariablesFromOldDist_preCollison. Preference for here?
6767 this->initRunCorrection();
6770 // set number outer Bnd cells for lblpt coupler
6771 m_noOuterBndryCells = this->m_bndCnd->m_bndCells.size();
6772 // post-propagation boundary conditions
6773 RECORD_TIMER_START(m_t.propagationBC);
6774 postPropagationBc();
6775 RECORD_TIMER_STOP(m_t.propagationBC);
6779 if(m_updateMacroscopicLocation == POSTPROPAGATION) {
6780 updateVariablesFromOldDist();
6783 RECORD_TIMER_STOP(m_t.solutionStep);
6789void LbSolver<nDim>::writeInfo() {
6792 if(!isActive()) return;
6793 RECORD_TIMER_START(m_t.residual);
6795 if(globalTimeStep % m_residualInterval == 0) {
6798 RECORD_TIMER_STOP(m_t.residual);
6803void LbSolver<nDim>::initSolver() {
6806 // set cell properties
6807 resetActiveCellList();
6809 // reset interface BC cell interpolation neighbors
6810 // TODO labels:LB this can not be called by treatInterfaceCell since a_isBndry is not
6811 // set correctly in advance.
6812 correctInterfaceBcCells();
6814 if(m_externalForcing) {
6815 initPressureForce();
6817 if(m_EELiquid.gravity || m_externalForcing) {
6821 if(grid().isActive()) {
6822 initSrcTermController();
6823 (this->*m_initializeMethod)();
6824 m_bndCnd->initializeBcData();
6828 m_log << endl << endl;
6839void LbSolver<nDim>::prepareAdaptation() {
6841 if(!m_refinedParents.empty()) m_refinedParents.clear();
6849void LbSolver<nDim>::setSensors(std::vector<std::vector<MFloat>>& sensors, std::vector<MFloat>& sensorWeight,
6850 std::vector<std::bitset<64>>& sensorCellFlag, std::vector<MInt>& sensorSolverId) {
6851 cerr0 << "Set
" << this->m_noSensors << " Sensors
for Adaptation
" << endl;
6854 ASSERT(m_freeIndices.empty(), "");
6856 // The sensors are added at the end of the previous sensor-vectors!
6857 // Implement sensor in the for loop below in this order.
6858 const auto sensorOffset = (signed)sensors.size();
6859 ASSERT(sensorOffset == 0 || grid().raw().treeb().noSolvers() > 1, "");
6860 sensors.resize(sensorOffset + this->m_noSensors, vector<MFloat>(grid().raw().m_noInternalCells, F0));
6861 sensorWeight.resize(sensorOffset + this->m_noSensors, -1);
6862 sensorCellFlag.resize(grid().raw().m_noInternalCells, sensorOffset + this->m_noSensors);
6863 sensorSolverId.resize(sensorOffset + this->m_noSensors, solverId());
6864 ASSERT(sensorOffset + this->m_noSensors < CartesianGrid<nDim>::m_maxNoSensors, "Increase bitset size!
");
6867 for(MInt sen = 0; sen < this->m_noSensors; sen++) {
6868 sensorWeight[sensorOffset + sen] = this->m_sensorWeight[sen];
6873 // only set sensors if the adaptation was globally triggered in buildLevelSetTubeCG!
6875 // Perform exchange to update Halos that are involved in gradient computation
6876 if(noNeighborDomains() > 0) {
6878 exchangeOldDistributions();
6881 if(domainId() == 0) {
6882 cerr << "Setting sensor
for the lb-solver adaptation!
" << endl;
6885 for(MInt sen = 0; sen < this->m_noSensors; sen++) {
6886 (this->*(this->m_sensorFnPtr[sen]))(sensors, sensorCellFlag, sensorWeight, sensorOffset, sen);
6897void LbSolver<nDim>::getSolverSamplingProperties(std::vector<MInt>& samplingVarIds,
6898 std::vector<MInt>& noSamplingVars,
6899 std::vector<std::vector<MString>>& samplingVarNames,
6900 const MString featureName) {
6903 // Read sampling variable names
6904 std::vector<MString> varNamesList;
6905 MInt noVars = readSolverSamplingVarNames(varNamesList, featureName);
6907 // Set default sampling variables if none specified
6909 varNamesList.emplace_back("LB_PV");
6913 for(MInt i = 0; i < noVars; i++) {
6914 const MInt samplingVar = string2enum(varNamesList[i]);
6915 std::vector<MString> varNames;
6917 auto samplingVarIt = std::find(samplingVarIds.begin(), samplingVarIds.end(), samplingVar);
6918 if(samplingVarIt != samplingVarIds.end()) {
6919 TERMM(1, "Sampling variable
'" + varNamesList[i] + "' already specified.
");
6922 switch(samplingVar) {
6924 ASSERT(m_noVariables == nDim + 1,
6925 "Error: additional primitive variables not supported
for sampling data feature yet
");
6927 samplingVarIds.push_back(LB_PV);
6928 noSamplingVars.push_back(m_noVariables + 1); //+1 because of addition of p
6930 varNames.resize(m_noVariables + 1); //+1 because of addition of p
6931 varNames[PV->U] = "u
";
6932 varNames[PV->V] = "v
";
6933 if constexpr(nDim == 3) {
6934 varNames[PV->W] = "w
";
6936 varNames[PV->P] = "p
"; // add a p variable to be sampled, that can be induced from rho
6937 varNames[PV->RHO] = "rho
";
6939 samplingVarNames.push_back(varNames);
6943 TERMM(1, "Unknown sampling variable:
" + varNamesList[i]);
6955void LbSolver<nDim>::finalizeAdaptation() {
6958 if(!grid().isActive()) return;
6959 if(m_maxNoSets > -1) {
6960 mDeallocate(m_sendBufferMB);
6961 mDeallocate(m_receiveBufferMB);
6962 mDeallocate(m_associatedBodyIds);
6963 mDeallocate(m_levelSetValues);
6964 mDeallocate(m_G0CellMapping);
6965 mDeallocate(m_isG0CandidateOfSet);
6967 // TODO: noNeighborDomains > 0 instead?
6968 if(noDomains() > 1) {
6971 for(MInt d = 0; d < noNeighborDomains(); d++) {
6972 sumwin += noWindowCells(d);
6973 sumhalo += noHaloCells(d);
6975 mAlloc(m_sendBufferMB, sumwin, "m_sendBufferMB
", 0, AT_);
6976 mAlloc(m_receiveBufferMB, sumhalo, "m_receiveBufferMB
", 0, AT_);
6980 mAlloc(m_associatedBodyIds, m_noLevelSetsUsedForMb * a_noCells(), "m_associatedBodyIds
", -1, AT_);
6981 mAlloc(m_levelSetValues, m_noLevelSetsUsedForMb * a_noCells(), "m_levelSetValues
", F0, AT_);
6982 mAlloc(m_G0CellMapping, a_noCells(), "m_G0CellMapping
", -1, AT_);
6984 MInt startSet = m_levelSetId; // if bodyToSetMode = 11 the 0th set is the collected set, which should be skiped
6985 MInt arrayLength = m_maxNoSets - startSet;
6986 mAlloc(m_isG0CandidateOfSet, a_noCells(), arrayLength, "m_isG0CandidateOfSet
", false, AT_);
6988 // TODO: noNeighborDomains > 0 instead?
6989 if(noDomains() > 1) {
6990 grid().updateLeafCellExchange();
6993 resetActiveCellList();
6996 // These status flag distinguishes if an adaptation was performed or not.
6997 // If adaptation was performed, recalcIds need to be passed to saveUVWRhoTOnlyPar
6998 // because the ids in the newly written grid file differ from the
6999 // ids used in solver/grid.
7000 m_adaptationSinceLastRestart = true;
7001 m_adaptationSinceLastSolution = true;
7017void LbSolver<nDim>::postAdaptation() {
7020 for(MInt i = 0; i < grid().raw().noNeighborDomains(); i++) {
7021 for(MInt j = 0; j < grid().raw().noHaloCells(i); j++) {
7022 MInt gridId = grid().raw().m_haloCells[i][j];
7023 MInt l_solverId = this->grid().tree().grid2solver(gridId);
7024 if(l_solverId < 0) {
7028 if(!grid().raw().treeb().solver(gridId, m_solverId)) {
7029 cerr << "removing
a Halo-cell
" << endl;
7030 this->removeCellId(l_solverId);
7035 // In meshAdaptation() the order of window and halo cells can change
7036 // even if no adaptation happens
7037 this->compactCells();
7039 if(!g_multiSolverGrid) {
7040 for(MInt gridCellId = 0; gridCellId < grid().raw().treeb().size(); gridCellId++) {
7041 ASSERT(grid().tree().solver2grid(gridCellId) == gridCellId, "");
7042 ASSERT(grid().tree().grid2solver(gridCellId) == gridCellId, "");
7046 grid().updateOther();
7047 updateDomainInfo(grid().domainId(), grid().noDomains(), grid().mpiComm(), AT_);
7049 if(!grid().isActive()) return;
7051 m_cells.size(c_noCells());
7052 m_freeIndices.clear();
7054 if(!g_multiSolverGrid) ASSERT(a_noCells() == c_noCells() && c_noCells() == grid().raw().treeb().size(), "");
7057 // set the cell weights.... i dont know why this is done by the solvers right now
7058 for(MInt i = 0; i < grid().raw().treeb().size(); i++) {
7059 grid().raw().treeb().weight(i) = 1;
7062 updateCellCollectorFromGrid();
7064 grid().findEqualLevelNeighborsParDiagonal(false);
7066 m_isRefined = (grid().maxUniformRefinementLevel() < maxLevel());
7071 prepareCommunicationReduced();
7073 prepareCommunicationNormal();
7078 // Do an exchange to have the right values in halo cells
7079 if(noNeighborDomains() > 0) {
7081 exchangeOldDistributions();
7088 if(!grid().isActive()) return;
7089 if(m_maxNoSets > -1) {
7090 mDeallocate(m_associatedBodyIds);
7091 mDeallocate(m_levelSetValues);
7092 mDeallocate(m_G0CellMapping);
7093 mDeallocate(m_isG0CandidateOfSet);
7095 mAlloc(m_associatedBodyIds, m_noLevelSetsUsedForMb * a_noCells(), "m_associatedBodyIds
", -1, AT_);
7096 mAlloc(m_levelSetValues, m_noLevelSetsUsedForMb * a_noCells(), "m_levelSetValues
", F0, AT_);
7097 mAlloc(m_G0CellMapping, a_noCells(), "m_G0CellMapping
", -1, AT_);
7099 MInt startSet = m_levelSetId; // if bodyToSetMode = 11 the 0th set is the collected set, which should be skiped
7100 MInt arrayLength = m_maxNoSets - startSet;
7101 mAlloc(m_isG0CandidateOfSet, a_noCells(), arrayLength, "m_isG0CandidateOfSet
", false, AT_);
7104 // Initialize refined cells after solver is updated because spatial interpolation
7105 // accross domain interfaces can become necessary. Also see refineCell.
7106 if(globalTimeStep < 1) {
7107 resetActiveCellList(1); // needs to be done here, since now bndry cells are known
7108 (this->*m_initializeMethod)();
7109 m_bndCnd->initializeBcData();
7111 initializeRefinedCellsPerLevel();
7112 resetActiveCellList(1);
7113 m_bndCnd->initializeBcData();
7114 initializeNewInterfaceParents();
7117 // New cells have been initialized, reset flag
7118 for(MInt i = 0; i < a_noCells(); i++) {
7119 a_hasProperty(i, SolverCell::WasNewlyCreated) = false;
7122 ASSERT(m_freeIndices.empty(), "");
7136void LbSolver<nDim>::refineCell(const MInt gridCellId) {
7137 // Get solver cell id of the parent
7138 MInt solverParentId = grid().tree().grid2solver(gridCellId);
7139 if(!g_multiSolverGrid) ASSERT(solverParentId == gridCellId, "");
7141 ASSERT(grid().raw().a_hasProperty(gridCellId, Cell::WasRefined), "");
7143 // Loop over newly created children of the parent cell
7144 for(MInt child = 0; child < grid().m_maxNoChilds; child++) {
7145 // Get grid cell id of new child
7146 MInt gridChildId = grid().raw().a_childId(gridCellId, child);
7147 if(gridChildId == -1) continue;
7149 if(!g_multiSolverGrid) ASSERT(grid().raw().a_hasProperty(gridChildId, Cell::WasNewlyCreated), "");
7151 // If solver is inactive all cells musst be halo cells!
7152 if(!grid().isActive()) ASSERT(grid().raw().a_isHalo(gridChildId), "");
7153 // If child exists in grid but is not located inside solver geometry
7154 if(!grid().solverFlag(gridChildId, solverId())) continue;
7156 // Create cell id for new child in solver
7157 const MInt solverChildId = this->createCellId(gridChildId);
7158 a_hasProperty(solverChildId, SolverCell::WasNewlyCreated) = true;
7160 m_refinedParents.insert(solverParentId);
7171void LbSolver<nDim>::removeChilds(const MInt gridCellId) {
7174 // Get solver id of the parent that will be coarsen
7175 MInt solverParentId = grid().tree().grid2solver(gridCellId);
7176 ASSERT(solverParentId > -1 && solverParentId < m_cells.size(), "solverParentId is:
" << solverParentId);
7177 ASSERT(c_noChildren(solverParentId) > 0, "");
7178 if(!g_multiSolverGrid) ASSERT(solverParentId == gridCellId, "");
7180 // SOLVER SPECIFIC PART
7181 // Initialize variables of parent from its children
7182 this->removeChildsLb(solverParentId);
7184 // Remove solver ids of deleted children
7185 for(MInt c = 0; c < grid().m_maxNoChilds; c++) {
7186 MInt childId = c_childId(solverParentId, c);
7187 if(childId < 0) continue;
7188 this->removeCellId(childId);
7190 if(!g_multiSolverGrid) {
7191 ASSERT((grid().raw().treeb().size() - m_cells.size()) <= grid().m_maxNoChilds, "");
7199void LbSolver<nDim>::removeCell(const MInt gridCellId) {
7202 // Get solver id of the parent that will be coarsen
7203 MInt cellId = grid().tree().grid2solver(gridCellId);
7205 ASSERT(cellId < m_cells.size(), "cellId is:
" << cellId);
7206 ASSERT(c_noChildren(cellId) == 0, "");
7208 this->removeCellId(cellId);
7212// this function should be moved to Solver as soon as cartesiansolver.h has been removed!!!
7213// this function should be moved to Solver as soon as cartesiansolver.h has been removed!!!
7214// this function should be moved to Solver as soon as cartesiansolver.h has been removed!!!
7216void LbSolver<nDim>::resizeGridMap() {
7217 grid().resizeGridMap(m_cells.size());
7226void LbSolver<nDim>::swapCells(const MInt cellId0, const MInt cellId1) {
7227 if(cellId1 == cellId0) return;
7229 for(MInt v = 0; v < m_noVariables; v++) {
7230 std::swap(a_variable(cellId1, v), a_variable(cellId0, v));
7231 std::swap(a_oldVariable(cellId1, v), a_oldVariable(cellId0, v));
7233 for(MInt dir = 0; dir < m_noDistributions; dir++) {
7234 std::swap(a_distribution(cellId1, dir), a_distribution(cellId0, dir));
7235 std::swap(a_oldDistribution(cellId1, dir), a_oldDistribution(cellId0, dir));
7238 for(MInt dir = 0; dir < m_noDistributions; dir++) {
7239 std::swap(a_distributionThermal(cellId1, dir), a_distributionThermal(cellId0, dir));
7240 std::swap(a_oldDistributionThermal(cellId1, dir), a_oldDistributionThermal(cellId0, dir));
7244 for(MInt dir = 0; dir < m_noDistributions; dir++) {
7245 std::swap(a_distributionTransport(cellId1, dir), a_distributionTransport(cellId0, dir));
7246 std::swap(a_oldDistributionTransport(cellId1, dir), a_oldDistributionTransport(cellId0, dir));
7249 std::swap(m_cells.allProperties(cellId1), m_cells.allProperties(cellId0));
7250 std::swap(a_nu(cellId1), a_nu(cellId0));
7251 std::swap(a_kappa(cellId1), a_kappa(cellId0));
7252 std::swap(a_diffusivity(cellId1), a_diffusivity(cellId0));
7254 std::swap(a_bndId(cellId1), a_bndId(cellId0));
7255 std::swap(a_level(cellId1), a_level(cellId0));
7257 if(!m_refinedParents.empty()) {
7258 auto it0 = m_refinedParents.find(cellId0);
7259 auto it1 = m_refinedParents.find(cellId1);
7260 if(it0 != m_refinedParents.end() && it1 == m_refinedParents.end()) {
7261 // nothing to be done
7262 } else if(it0 != m_refinedParents.end()) {
7263 m_refinedParents.erase(it0);
7264 m_refinedParents.insert(cellId1);
7265 } else if(it1 != m_refinedParents.end()) {
7266 m_refinedParents.erase(it1);
7267 m_refinedParents.insert(cellId0);
7278void LbSolver<nDim>::swapProxy(const MInt cellId0, const MInt cellId1) {
7279 grid().swapGridIds(cellId0, cellId1);
7288void LbSolver<nDim>::resetComm() {
7290 mDeallocate(m_noWindowDistDataPerDomain);
7291 mDeallocate(m_noHaloDistDataPerDomain);
7292 mDeallocate(m_nghbrOffsetsWindow);
7293 mDeallocate(m_nghbrOffsetsHalo);
7294 mDeallocate(m_sendBuffers);
7295 mDeallocate(m_receiveBuffers);
7296 mDeallocate(m_needsFurtherExchange);
7297 mDeallocate(m_windowDistsForExchange);
7298 mDeallocate(m_haloDistsForExchange);
7299 if(!m_nonBlockingComm) {
7300 mDeallocate(mpi_request);
7302 mDeallocate(mpi_requestS);
7303 mDeallocate(mpi_requestR);
7306 mDeallocate(m_dataBlockSizes);
7307 mDeallocate(m_baseAddresses);
7308 mDeallocate(m_nghbrOffsetsWindow);
7309 mDeallocate(m_nghbrOffsetsHalo);
7310 mDeallocate(m_sendBuffers);
7311 mDeallocate(m_receiveBuffers);
7312 if(!m_nonBlockingComm) {
7313 mDeallocate(mpi_request);
7315 mDeallocate(mpi_requestS);
7316 mDeallocate(mpi_requestR);
7326void LbSolver<nDim>::resetCellLists(MBool resize /*= true*/) {
7328 // resize list of active cells
7329 mDeallocate(m_activeCellList);
7330 mAlloc(m_activeCellList, grid().noCells(), "m_activeCellList
", AT_);
7334 resetInterfaceCells();
7335 initializeInterfaceCells();
7336 buildInterfaceCells();
7337 setInterpolationNeighbors();
7338 setInterpolationCoefficients();
7348MBool LbSolver<nDim>::prepareRestart(MBool writeRestart, MBool& writeGridRestart) {
7351 writeGridRestart = false;
7353 if(((globalTimeStep % m_restartInterval) == 0) || writeRestart) {
7354 writeRestart = true;
7356 if(m_adaptationSinceLastRestart) {
7357 writeGridRestart = true;
7361 return writeRestart;
7370void LbSolver<nDim>::reIntAfterRestart(MBool doneRestart) {
7374 m_adaptationSinceLastRestart = false;
7383void LbSolver<nDim>::writeRestartFile(const MBool writeRestart, const MBool /*writeBackup*/, const MString gridFileName,
7384 MInt* recalcIdTree) {
7388 stringstream fileName;
7389 fileName << outputDir() << "restart_
" << getIdentifier() << globalTimeStep << ParallelIo::fileExt();
7391 std::vector<MInt> recalcCellIdsSolver(0);
7393 MInt noInternalCellIds;
7394 std::vector<MInt> reorderedCellIds(0);
7395 this->calcRecalcCellIdsSolver(recalcIdTree, noCells, noInternalCellIds, recalcCellIdsSolver, reorderedCellIds);
7396 MInt* pointerRecalcIds = (recalcIdTree == nullptr) ? nullptr : recalcCellIdsSolver.data();
7397 saveRestartWithDistributionsPar(fileName.str().c_str(), gridFileName.c_str(), pointerRecalcIds);
7406void LbSolver<nDim>::setCellWeights(MFloat* solverCellWeight) {
7408 const MInt noCellsGrid = grid().raw().treeb().size();
7409 const MInt offset = noCellsGrid * solverId();
7411 for(MInt cellId = 0; cellId < a_noCells(); cellId++) {
7412 const MInt gridCellId = grid().tree().solver2grid(cellId);
7413 const MInt id = gridCellId + offset;
7414 solverCellWeight[id] = F1; // 0.2;
7419void LbSolver<nDim>::initializeNewInterfaceParents() {
7420 const MInt interfaceLevels = maxLevel() - minLevel();
7422 for(MInt level = 0; level < interfaceLevels; level++) {
7423 for(MInt p = 0; p < m_noInterfaceParents[level]; p++) {
7424 auto& parent = m_interfaceParents[level]->a[p];
7425 const MInt parentId = parent.m_cellId;
7426 if(!a_wasActive(parentId) && a_isActive(parentId)) {
7427 this->removeChildsLb(parentId);
7443void LbSolver<nDim>::initializeRefinedCellsPerLevel() {
7444 for(MInt level = grid().maxUniformRefinementLevel(); level < grid().maxLevel(); level++) {
7445 // Exchange to update Halo cells that are needed for spatial interpolation
7446 // TODO labels:LB Find out how to exchange cells only on the current level!
7447 if(noNeighborDomains() > 0) {
7449 exchangeOldDistributions();
7452 for(auto&& parentId : m_refinedParents) {
7453 if(c_level(parentId) != level) continue;
7454 if(a_isHalo(parentId)) continue;
7455 if(!g_multiSolverGrid)
7456 if(!grid().raw().a_hasProperty(parentId, Cell::WasRefined)) continue;
7458 // Gather children of the refined parent cell and call initialization method
7459 std::vector<MInt> childIds(grid().m_maxNoChilds, -1);
7460 for(MInt child = 0; child < grid().m_maxNoChilds; child++) {
7461 MInt childId = c_childId(parentId, child);
7465 if(!g_multiSolverGrid) ASSERT(grid().raw().a_hasProperty(childId, Cell::WasNewlyCreated), "");
7467 childIds[child] = childId;
7469 this->refineCellLb(parentId, childIds.data());
7475void LbSolver<nDim>::initSolverSamplingVariables(const std::vector<MInt>& varIds,
7476 const std::vector<MInt>& noSamplingVars) {
7479 // TODO labels:LB,DLB enable use with balancing
7480 if(m_isInitSamplingVars) {
7481 m_log << "Sampling variables already initialized.
" << std::endl;
7485 for(MUint i = 0; i < varIds.size(); i++) {
7486 MFloat** varPointer = nullptr;
7487 MInt dataBlockSize = noSamplingVars[i];
7489 // No additional storage for primitive variables
7490 if(varIds[i] == LB_PV) {
7494 if(dataBlockSize > 0) {
7495 mAlloc(varPointer, a_noCells(), dataBlockSize, "m_samplingVariables_
" + std::to_string(varIds[i]), 0.0, AT_);
7497 m_samplingVariables.push_back(varPointer);
7500 // Set flag to avoid multiple initializations when using multiple sampling features
7501 m_isInitSamplingVars = true;
7505void LbSolver<nDim>::calcSamplingVariables(const std::vector<MInt>& varIds, const MBool NotUsed(exchange)) {
7506 for(MUint i = 0; i < varIds.size(); i++) {
7513 TERMM(1, "Invalid variable
id");
7521void LbSolver<nDim>::interpolateVariablesInCell(const MInt cellId, const MFloat* /*position*/,
7522 MFloat* interpolationResult) {
7523 // TODO labels:LB Add trilinear or least square interpolation
7524 interpolationResult[PV->RHO] = a_variable(cellId, PV->RHO);
7525 for(MInt d = 0; d < nDim; d++) {
7526 interpolationResult[PV->U + d] = a_variable(cellId, PV->U + d);
7528 interpolationResult[PV->P] = interpolationResult[PV->RHO] * CSsq;
7533void LbSolver<nDim>::calcSamplingVarAtPoint(const MFloat* point, const MInt id, const MInt sampleVarId, MFloat* state,
7534 const MBool NotUsed(interpolate)) {
7535 // Note: interpolateVariablesInCell() requires variables for all cells
7536 switch(sampleVarId) {
7538 interpolateVariablesInCell(id, point, state);
7542 TERMM(1, "Invalid variable
id");
7550void LbSolver<nDim>::saveSolverSolution(MBool /*forceOutput*/, const MBool /*finalTimeStep*/) {
7551 if(globalTimeStep % m_solutionInterval == 0 && globalTimeStep >= m_solutionOffset) saveOutput();
7553 if(m_fftInterval > 0 && globalTimeStep % m_fftInterval == 0
7555 RECORD_TIMER_START(m_t.fft);
7556 computeFFTStatistics();
7557 RECORD_TIMER_STOP(m_t.fft);
7562void LbSolver<nDim>::finalizeInitSolver() {
7563 if(!isActive()) return;
7564 if(m_updateMacroscopicLocation == POSTPROPAGATION && !m_restartFile) {
7565 updateVariablesFromOldDist();
7568 outputInitSummary();
7577void LbSolver<nDim>::outputInitSummary() {
7578 using namespace maia::logtable;
7580 // TODO labels:LB,IO This should be moved into an appropriate sysEqn-like structure...
7581 const MInt noVars = nDim + 1;
7582 std::vector<MString> cvs{};
7583 cvs.push_back("rho
");
7584 cvs.push_back("rho_u
");
7585 cvs.push_back("rho_v
");
7587 cvs.push_back("rho_w
");
7590 // INIT SUMMARY FRAME
7591 Frame summary(" SOLVER
" + std::to_string(solverId()) + " INITIALIZATION SUMMARY AT TIME STEP
"
7592 + std::to_string(globalTimeStep));
7594 // PROBLEM SUMMARY GROUP
7595 Group& problem = summary.addGroup(" PROBLEM SUMMARY
");
7596 problem.addData("System of Equations
", solverMethod());
7597 problem.addData("Number of dimensions
", nDim);
7598 Data& vars = problem.addData("Number of variables
", noVars);
7601 for(auto&& cv : cvs) {
7603 vars.addData("Conservative variable name(s)
", cv);
7605 vars.addData("", cv);
7610 problem.addData("Reynolds Number
", m_Re);
7611 problem.addData("Mach Number
", m_Ma);
7612 problem.addData("Reference length (
LB)
", m_referenceLength);
7614 problem.addData("Restart
", m_restartFile);
7616 const MString restartFileName =
7617 restartDir() + "restart_
" + getIdentifier() + std::to_string(m_restartTimeStep) + ParallelIo::fileExt();
7618 problem.addData("Initial condition
", restartFileName);
7620 Data& initData = problem.addData("Initial condition
", m_initMethod);
7622 initData.addData("tanhInit Re
", m_initRe);
7623 initData.addData("tanhInit start time
", m_initStartTime);
7624 initData.addData("tanhInit time
", m_initTime);
7627 problem.addData("Start time (non-dimensionalized)
", globalTimeStep);
7628 problem.addData("Final time (non-dimensionalized)
", globalTimeStep + g_timeSteps);
7630 // DISCRETIZATION SUMMARY GROUP
7631 Group& discret = summary.addGroup("DISCRETIZATION SUMMARY
");
7632 discret.addData("Number of distributions
", m_noDistributions);
7633 discret.addData("Collision
operator", solverMethod());
7635 Data& omega = discret.addData("Collision frequency
", m_omega);
7636 omega.addData("Relaxation time
", 1.0 / m_omega);
7637 omega.addData("Num. viscosity
", m_nu);
7639 Data& omegaT = discret.addData("Collision frequency thermal
", m_omegaT);
7640 omegaT.addData("Relaxation time thermal
", 1.0 / m_omegaT);
7641 omegaT.addData("Num. diffusion thermal
", m_kappa);
7644 Data& omegaD = discret.addData("Collision frequency transport
", m_omegaD);
7645 omegaD.addData("Relaxation time transport
", 1.0 / m_omegaD);
7646 omegaD.addData("Num. diffusion transport
", m_diffusivity);
7649 discret.addData("External forcing
", m_externalForcing);
7653 interfaceMethod = Context::getSolverProperty<MString>("interfaceMethod
", m_solverId, AT_, &interfaceMethod);
7655 discret.addData("Interface method
", interfaceMethod);
7657 discret.addData("Adaptation init method
", m_adaptationInitMethod);
7659 // PARALLELIZATION SUMMARY
7660 Group& parallel = summary.addGroup("PARALLELIZATION SUMMARY
");
7661 parallel.addData("Domain
id", domainId());
7662 parallel.addData("Number of neighbor domains
", grid().noNeighborDomains());
7663 parallel.addData("Total number of domains
", noDomains());
7665 // GRID SUMMARY LOCAL
7666 Group& gridLocal = summary.addGroup("GRID SUMMARY (
LOCAL)
");
7667 Data& minlvl = gridLocal.addData("Minimum grid level
", grid().minLevel());
7668 minlvl.addData("cell length
", c_cellLengthAtLevel(minLevel()));
7670 Data& maxlvl = gridLocal.addData("Maximum grid level
", grid().maxLevel());
7671 maxlvl.addData("cell length
", c_cellLengthAtLevel(maxLevel()));
7673 gridLocal.addBlank();
7674 Data& noCellsLocal = gridLocal.addData("Number of cells
", a_noCells());
7675 noCellsLocal.addData("internal cells
", grid().noInternalCells());
7676 noCellsLocal.addData("halo cells
", a_noCells() - grid().noInternalCells());
7677 gridLocal.addData("Number of active cells
", m_currentMaxNoCells);
7678 gridLocal.addData("Max. number of cells (collector size)
", grid().raw().treeb().capacity());
7679 gridLocal.addData("Memory utilization
", 100.0 * grid().noCells() / grid().raw().treeb().capacity());
7680 gridLocal.addBlank();
7682 const MInt levelRange = maxLevel() - minLevel();
7683 const MInt totalNoInterfaceChildren =
7684 std::accumulate(&m_noInterfaceChildren[0], &m_noInterfaceChildren[levelRange], 0);
7685 Data& ic = gridLocal.addData("Number of interface children
", totalNoInterfaceChildren);
7686 for(MInt i = 0; i < levelRange; i++) {
7687 const MInt noInterfaceChildren = m_noInterfaceChildren[i];
7688 if(noInterfaceChildren > 0) {
7689 ic.addData("on level
" + std::to_string(minLevel() + i + 2), noInterfaceChildren);
7692 const MInt totalNoInterfaceParents =
7693 std::accumulate(&m_noInterfaceParents[0], &m_noInterfaceParents[levelRange], 0);
7694 Data& ip = gridLocal.addData("Number of interface parents
", totalNoInterfaceParents);
7695 for(MInt i = 0; i < levelRange; i++) {
7696 const MInt noInterfaceParents = m_noInterfaceParents[i];
7697 if(noInterfaceParents > 0) {
7698 ip.addData("on level
" + std::to_string(minLevel() + i + 1), noInterfaceParents);
7703 // Finally build string for summary frame
7704 MString s = summary.buildString();
7706 if(domainId() == 0) {
7709 m_log << s << std::endl;
7714void LbSolver<nDim>::preTimeStep() {
7717 // update m_time for output
7722//#################################################################################
7723// MOVING BOUNDARY PART STARTS HERE!!!!!!
7724//#################################################################################
7727void LbSolver<nDim>::preCoupleLs(std::vector<MInt>& maxGCellLevels) {
7730 if(!grid().isActive()) {
7734 RECORD_TIMER_START(m_t.findG0Cells);
7736 this->exchangeData(&a_levelSetFunctionMB(0, 0));
7738 MInt startSet = m_levelSetId; // if bodyToSetMode = 11 the 0th set is the collected set, which should be skiped
7740 resetActiveCellList(2);
7744 prepareCommunicationReduced();
7747 RECORD_TIMER_START(m_t.findG0Candidates);
7748 findG0Candidates(maxGCellLevels); // findes all the G0Candidates of all Sets
7749 RECORD_TIMER_STOP(m_t.findG0Candidates);
7751 MBool gapClosure = false;
7752 m_G0Candidates.clear();
7754 MInt candidates = 0;
7755 for(MInt i = 0; i < a_noCells(); i++) {
7756 MBool isG0Candidate = false;
7757 for(MInt set = 0; set < (m_maxNoSets - startSet); set++) {
7758 if(a_isG0CandidateOfSet(i, set)) {
7759 isG0Candidate = true;
7763 m_G0CellMapping[i] = -1;
7765 if(a_isHalo(i)) continue;
7766 m_G0Candidates.emplace_back();
7767 m_G0Candidates[candidates].cellId = i;
7768 m_G0CellMapping[i] = candidates;
7773 MBoolScratchSpace isGapCell(a_noCells(), AT_, "isGapCell
");
7774 isGapCell.fill(false);
7776 RECORD_TIMER_START(m_t.geomNodal);
7777 m_geometryIntersection->computeNodalValues(m_G0Candidates, &m_G0CellMapping[0], &a_levelSetFunctionMB(0, 0),
7778 &a_associatedBodyIds(0, 0), &isGapCell(0), gapClosure);
7779 RECORD_TIMER_STOP(m_t.geomNodal);
7781 if(noNeighborDomains() > 0) {
7782 RECORD_TIMER_START(m_t.geomExchange);
7784 // Save pointer to leafWindow/leafHaloCells per neighborDomain to get 'const MInt**' argument
7785 ScratchSpace<const MInt*> p_leafWindowCells(noNeighborDomains(), AT_, "p_leafWindowCells
");
7786 ScratchSpace<const MInt*> p_leafHaloCells(noNeighborDomains(), AT_, "p_leafHaloCells
");
7787 ScratchSpace<MInt> noLeafWindowCells(noNeighborDomains(), AT_, "noLeafWindowCells
");
7788 ScratchSpace<MInt> noLeafHaloCells(noNeighborDomains(), AT_, "noLeafHaloCells
");
7789 for(MInt d = 0; d < noNeighborDomains(); d++) {
7790 p_leafWindowCells[d] = &grid().leafWindowCell(d, 0);
7791 p_leafHaloCells[d] = &grid().leafHaloCell(d, 0);
7792 noLeafWindowCells[d] = grid().noLeafWindowCells(d);
7793 noLeafHaloCells[d] = grid().noLeafHaloCells(d);
7796 m_geometryIntersection->exchangeNodalValues(p_leafWindowCells.data(), noLeafWindowCells.data(),
7797 p_leafHaloCells.data(), m_G0Candidates, &m_G0CellMapping[0]);
7798 RECORD_TIMER_STOP(m_t.geomExchange);
7801 if(m_G0CellList != nullptr) {
7802 mDeallocate(m_G0CellList);
7804 if(m_nodalGValues != nullptr) {
7805 mDeallocate(m_nodalGValues);
7808 if(!m_G0Candidates.empty()) {
7809 MInt noNodalGValues = m_noLevelSetsUsedForMb * (IPOW3[nDim] - 1);
7810 m_noG0CandidatesTotal = (MInt)m_G0Candidates.size();
7811 mAlloc(m_G0CellList, m_noG0CandidatesTotal, "m_G0CellList
", -1, AT_);
7812 mAlloc(m_nodalGValues, m_noG0CandidatesTotal, noNodalGValues, "m_nodalGValues
", F0, AT_);
7814 m_noG0CandidatesTotal = 0;
7818 std::fill_n(m_G0CellMapping, a_noCells(), -1);
7820 RECORD_TIMER_START(m_t.calcNodalValues);
7821 calcNodalLsValues();
7822 RECORD_TIMER_STOP(m_t.calcNodalValues);
7824 RECORD_TIMER_STOP(m_t.findG0Cells);
7828void LbSolver<nDim>::createBndryToBodyMapping(maia::coupling::Mapping& bndryToBodyMapping,
7829 maia::coupling::Mapping& bodyToBndryMapping) {
7830 // TODO labels:LB Use collector!
7831 // Create Mapping bndryCell -> bodyId
7832 bndryToBodyMapping.clear();
7833 bodyToBndryMapping.clear();
7834 for(MInt i = 0; i < m_currentNoG0Cells; i++) {
7835 const MInt bndryCellId = m_G0CellList[i];
7836 if(a_associatedBodyIds(bndryCellId, 0) >= 0) {
7837 bndryToBodyMapping.insert(i) = a_associatedBodyIds(bndryCellId, 0);
7838 bodyToBndryMapping.insert(a_associatedBodyIds(bndryCellId, 0)) = i;
7843 for(MInt candidate = 0; candidate < m_currentNoG0Cells; candidate++) {
7844 const MInt cellId = m_G0CellList[candidate];
7845 ASSERT(cellId > -1, "G0CellList broken!
");
7846 ASSERT(m_G0CellMapping[cellId] > -1, "G0CellMapping broken!
");
7847 std::stringstream err;
7848 err << "G0CellMapping broken
" << m_G0CellMapping[cellId] << " not eq to
" << candidate;
7849 ASSERT(m_G0CellMapping[cellId] == candidate, err.str());
7859void LbSolver<nDim>::resetActiveCellList(MInt mode) {
7861 //--store wasActive state-----------------------------------------------------
7862 for(MInt i = 0; i < a_noCells(); i++) {
7863 a_wasActive(i) = a_isActive(i);
7865 //--initialization of a_isActive----------------------------------------------
7867 setIsActiveDefaultStates();
7869 setInActiveBndryCells();
7871 setInActiveMBCells();
7873 // TODO: Between the above and the below part, coupler might be called in
7874 // future. Then the upper and lower part has to be split into several
7875 // functions. (resetActiveStates and resetActiveCellList)
7876 //--reset activeCellList------------------------------------------------------
7877 fillActiveCellList();
7880// Reset external forces vector for each cell
7882void LbSolver<nDim>::resetExternalSources() {
7884 maia::parallelFor(0, a_noCells(), [&](MInt cellId) {
7885 for(MInt dim = 0; dim < nDim; dim++) {
7886 a_externalForces(cellId, dim) = 0.0;
7892void LbSolver<nDim>::exchangeExternalSources() {
7895 ASSERT(m_particleMomentumCoupling, "");
7897 if(noDomains() > 1 && m_particleMomentumCoupling) {
7898 maia::mpi::reverseExchangeAddData(grid().neighborDomains(), grid().haloCells(), grid().windowCells(), mpiComm(),
7899 &a_externalForces(0, 0), nDim);
7909inline void LbSolver<nDim>::setIsActiveDefaultStates() {
7911 const MInt noCells = m_cells.size();
7913 for(MInt i = 0; i < noCells; i++) {
7914 const MBool l_isActive = c_isLeafCell(i) || a_isInterfaceParent(i);
7915 a_isActive(i) = (l_isActive) ? 1 : 0;
7918 for(MInt i = 0; i < noCells; i++) {
7919 const MBool l_isActive = c_isLeafCell(i);
7920 a_isActive(i) = (l_isActive) ? 1 : 0;
7923 if(m_solidLayerExtension) {
7924 if(m_initialActiveCells != nullptr) {
7925 mDeallocate(m_initialActiveCells);
7927 mAlloc(m_initialActiveCells, noCells, "m_initialActiveCells
", 1, AT_);
7928 for(MInt i = 0; i < noCells; i++) {
7930 a_isActive(i) = !m_geometry->pointIsInside2(&a_coordinate(i, 0));
7932 m_initialActiveCells[i] = a_isActive(i);
7942inline void LbSolver<nDim>::setInActiveBndryCells() {
7944 // find ids of wall boundaries
7945 vector<MInt> wallIds;
7946 for(MInt i = 0; i < (MInt)(m_bndCnd->m_bndCndSegIds.size()); i++) {
7947 MBool is_periodic = false;
7948 if(m_bndCnd->m_noPeriodicSegments != 0)
7949 for(MInt j = 0; j < m_bndCnd->m_noPeriodicSegments; j++)
7950 if(m_bndCnd->m_bndCndSegIds[i] == m_bndCnd->m_periodicSegmentsIds[j]) {
7954 MBool is_inout = false;
7955 for(MInt j = 0; j < m_bndCnd->m_noInOutSegments; j++)
7956 if(m_bndCnd->m_bndCndSegIds[i] == m_bndCnd->m_inOutSegmentsIds[j]) {
7960 if(!is_periodic & !is_inout) wallIds.push_back(i);
7962 // set non-fluid cell to inactive
7963 for(auto wallId : wallIds) {
7964 for(MInt i = m_bndCnd->m_bndCndOffsets[wallId]; i < m_bndCnd->m_bndCndOffsets[wallId + 1]; i++) {
7965 if(m_bndCnd->m_bndCells[i].m_isFluid == false) {
7966 const MInt id = m_bndCnd->m_bndCells[i].m_cellId;
7967 a_isActive(id) = false;
7979inline void LbSolver<nDim>::setInActiveMBCells() {
7981 if(m_maxNoSets == -1 || m_levelSetValues == nullptr) return;
7982 for(MInt i = 0; i < a_noCells(); i++) {
7983 if(m_solidLayerExtension) {
7984 if(!m_initialActiveCells[i]) continue;
7986 for(MInt set = 0; set < m_maxNoSets; set++) {
7987 if(a_associatedBodyIds(i, set) < 0) continue;
7988 if(!c_isLeafCell(i)) continue;
7989 if(a_levelSetFunctionMB(i, set) < 0) {
8003inline void LbSolver<nDim>::fillActiveCellList() {
8004 if(!isActive()) return;
8005 const MInt noLevels = maxLevel() - minLevel() + 1;
8007 // counting active cells - total and per level
8008 m_currentMaxNoCells = 0;
8009 ScratchSpace<MInt> noActiveCellsPerLevel(noLevels, AT_, "noActiveCellsPerLevel
");
8010 noActiveCellsPerLevel.fill(0);
8011 for(MInt i = 0; i < a_noCells(); i++) {
8013 const MInt lvl = a_level(i) - minLevel();
8014 m_activeCellList[m_currentMaxNoCells] = i;
8015 noActiveCellsPerLevel[lvl]++;
8016 m_currentMaxNoCells++;
8019 // fill left entries with non-valid value
8020 for(MInt i = m_currentMaxNoCells; i < a_noCells(); i++) {
8021 m_activeCellList[i] = -1;
8023 // m_log << " + number of active cells:
" << m_currentMaxNoCells << endl << endl;
8024 // fill offset for better accessing later
8025 if(m_activeCellListLvlOffset != nullptr) {
8026 mDeallocate(m_activeCellListLvlOffset);
8028 mAlloc(m_activeCellListLvlOffset, noLevels + 1, "m_activeCellListLvlOffset
", 0, AT_);
8029 m_activeCellListLvlOffset[0] = 0;
8030 for(MInt i = 1; i < noLevels + 1; i++) {
8031 m_activeCellListLvlOffset[i] = m_activeCellListLvlOffset[i - 1] + noActiveCellsPerLevel[i - 1];
8033 // sort active cells per level TODO labels:LB,toenhance sorting algorithm might be worth improving (HEAP-sort?)
8035 MInt level = minLevel();
8037 for(MInt l = 0; l < noLevels - 1; l++) {
8039 for(; i < m_activeCellListLvlOffset[l + 1]; i++) {
8040 if(a_level(m_activeCellList[i]) == level) {
8042 continue; // already at correct position
8044 for(; j < m_currentMaxNoCells; j++) {
8045 if(a_level(m_activeCellList[j]) == level) {
8046 std::swap(m_activeCellList[i], m_activeCellList[j]);
8047 break; // now with correct level, next position i
8063void LbSolver<nDim>::findG0Candidates(std::vector<MInt>& maxGCellLevels) {
8066 const MInt startSet = m_levelSetId;
8068 for(MInt set = startSet; set < m_maxNoSets; set++) {
8069 const MInt positionInArray = set - startSet;
8070 const MFloat maxDistanceG0 = 4.0 * c_cellLengthAtLevel(maxGCellLevels[set]);
8072 for(MInt i = 0; i < a_noCells(); i++) {
8073 a_isG0CandidateOfSet(i, positionInArray) = false;
8074 if(a_isHalo(i)) continue;
8075 if(!c_isLeafCell(i)) continue;
8076 if(a_associatedBodyIds(i, set) == -2) continue;
8077 // if a_associatedBodyId is set to -2 it means that the corresponding cell is
8078 // not cell belonging to the LS solver. Hence, no valid LS values and body Ids
8079 // can be given here. Cells which are located next to cells which are not
8080 // owned by the LS solver require special treatment. These cells should not
8081 // be G0 candidates because their nodel G-value cannot be calculated accurately.
8082 // Hence, these cells are filltered out.
8083 MBool hasNeighborInPureLbDomain = false;
8084 for(MInt n = 0; n < IPOW3[nDim] - 1; n++) {
8085 if(c_neighborId(i, n) < 0) {
8088 if(a_associatedBodyIds(c_neighborId(i, n), set) == -2) {
8089 hasNeighborInPureLbDomain = true;
8090 a_associatedBodyIds(i, set) = -1;
8095 if(hasNeighborInPureLbDomain) {
8096 a_isG0CandidateOfSet(i, positionInArray) = false;
8098 if(m_allowBndryAsG0) {
8099 if((fabs(a_levelSetFunctionMB(i, set)) < maxDistanceG0)) {
8100 a_isG0CandidateOfSet(i, positionInArray) = true;
8102 a_isG0CandidateOfSet(i, positionInArray) = false;
8105 if((fabs(a_levelSetFunctionMB(i, set)) < maxDistanceG0) && !a_isBndryCell(i)) {
8106 a_isG0CandidateOfSet(i, positionInArray) = true;
8108 a_isG0CandidateOfSet(i, positionInArray) = false;
8115 if(noDomains() > 1) {
8116 exchangeG0Candidates();
8128void LbSolver<nDim>::exchangeG0Candidates() {
8131 MInt startSet = m_levelSetId; // if bodyToSetMode = 11 the 0th set is the collected set, which should be skiped
8132 for(MInt set = 0; set < (m_maxNoSets - startSet); set++) {
8133 // 1. Gathered information (to which level set does the cell belong) of the window cells.
8134 MInt sendNeighOffset = 0;
8135 for(MInt n = 0; n < noNeighborDomains(); n++) {
8136 for(MInt j = 0; j < noWindowCells(n); j++) {
8137 m_sendBufferMB[sendNeighOffset + j] = a_isG0CandidateOfSet(windowCell(n, j), set);
8139 sendNeighOffset += noWindowCells(n);
8142 // 2. Send the gathered information.
8143 sendNeighOffset = 0;
8144 for(MInt n = 0; n < noNeighborDomains(); n++) {
8145 MInt bufsize = noWindowCells(n);
8146 MPI_Issend(&m_sendBufferMB[sendNeighOffset], bufsize, MPI_CXX_BOOL, neighborDomain(n), 0, mpiComm(),
8147 &mpi_request[n], AT_, "m_sendBufferMB[sendNeighOffset]
");
8148 sendNeighOffset += noWindowCells(n);
8151 // 3. Receive data from neighbors.
8153 MInt recNeighOffset = 0;
8154 for(MInt n = 0; n < noNeighborDomains(); n++) {
8155 MInt bufsize = noHaloCells(n);
8156 MPI_Recv(&m_receiveBufferMB[recNeighOffset], bufsize, MPI_CXX_BOOL, neighborDomain(n), 0, mpiComm(), &status, AT_,
8157 "m_receiveBufferMB[recNeighOffset]
");
8158 recNeighOffset += noHaloCells(n);
8161 for(MInt n = 0; n < noNeighborDomains(); n++)
8162 MPI_Wait(&mpi_request[n], &status, AT_);
8164 // 4. scatter the received data
8166 for(MInt n = 0; n < noNeighborDomains(); n++) {
8167 for(MInt j = 0; j < noHaloCells(n); j++) {
8168 a_isG0CandidateOfSet(haloCell(n, j), set) = m_receiveBufferMB[recNeighOffset + j];
8170 recNeighOffset += noHaloCells(n);
8184MInt LbSolver<nDim>::setUpLbInterpolationStencil(const MInt cellId, MInt* interpolationCells, MFloat* point) {
8187 std::function<MBool(const MInt, const MInt)> alwaysTrue = [&](const MInt, const MInt) { return true; };
8189 return this->setUpInterpolationStencil(cellId, interpolationCells, point, alwaysTrue, false);
8193MFloat LbSolver<nDim>::interpolateFieldDataLb(MInt* interpolationCells, MFloat* point, MInt set,
8194 std::function<MFloat(MInt, MInt)> scalarField,
8195 std::function<MFloat(MInt, MInt)> coordinate) {
8198 return this->template interpolateFieldData<true>(interpolationCells, point, set, scalarField, coordinate);
8202void LbSolver<nDim>::initializeMovingBoundaries() {
8205 //########################################################################
8206 // start moving boundary properties
8207 //########################################################################
8222 m_trackMovingBndry = true;
8223 m_trackMovingBndry = Context::getSolverProperty<MBool>("trackMovingBndry
", m_solverId, AT_, &m_trackMovingBndry);
8238 m_trackMbStart = Context::getSolverProperty<MInt>("trackMbStart
", m_solverId, AT_, &m_trackMbStart);
8252 m_trackMbEnd = numeric_limits<MInt>::max();
8253 m_trackMbEnd = Context::getSolverProperty<MInt>("trackMbEnd
", m_solverId, AT_, &m_trackMbEnd);
8255 if(m_sendBufferMB != nullptr) mDeallocate(m_sendBufferMB);
8256 if(m_receiveBufferMB != nullptr) mDeallocate(m_receiveBufferMB);
8257 if(noDomains() > 1 && grid().isActive()) {
8260 for(MInt d = 0; d < noNeighborDomains(); d++) {
8261 sumwin += noWindowCells(d);
8262 sumhalo += noHaloCells(d);
8264 mAlloc(m_sendBufferMB, sumwin, "m_sendBufferMB
", 0, AT_);
8265 mAlloc(m_receiveBufferMB, sumhalo, "m_receiveBufferMB
", 0, AT_);
8267 //########################################################################
8268 // end moving boundary properties
8269 //########################################################################
8271 //########################################################################
8272 // start updating ActiveCellList for Moving Boundaries
8273 //########################################################################
8276 if(m_associatedBodyIds != nullptr) mDeallocate(m_associatedBodyIds);
8277 mAlloc(m_associatedBodyIds, m_noLevelSetsUsedForMb * a_noCells(), "m_associatedBodyIds
", -1, AT_);
8278 if(m_levelSetValues != nullptr) mDeallocate(m_levelSetValues);
8279 mAlloc(m_levelSetValues, m_noLevelSetsUsedForMb * a_noCells(), "m_levelSetValues
", F0, AT_);
8281 MInt startSet = m_levelSetId; // if bodyToSetMode = 11 the 0th set is the collected set, which should be skiped
8282 MInt arrayLength = m_maxNoSets - startSet;
8283 if(m_isG0CandidateOfSet != nullptr) mDeallocate(m_isG0CandidateOfSet);
8284 mAlloc(m_isG0CandidateOfSet, a_noCells(), arrayLength, "m_isG0CandidateOfSet
", false, AT_);
8286 //########################################################################
8287 // end updating ActiveCellList for Moving Boundaries
8288 //########################################################################
8290 if(m_G0CellMapping != nullptr) mDeallocate(m_G0CellMapping);
8291 mAlloc(m_G0CellMapping, a_noCells(), "m_G0CellMapping
", -1, AT_);
8292 if(m_innerBandWidth != nullptr) mDeallocate(m_innerBandWidth);
8293 mAlloc(m_innerBandWidth, grid().maxRefinementLevel(), "m_innerBandWidth
", F0, AT_);
8294 if(m_bandWidth != nullptr) mDeallocate(m_bandWidth);
8295 mAlloc(m_bandWidth, grid().maxRefinementLevel(), "m_bandWidth
", 0, AT_);
8296 if(m_outerBandWidth != nullptr) mDeallocate(m_outerBandWidth);
8297 mAlloc(m_outerBandWidth, grid().maxRefinementLevel(), "m_outerBandWidth
", F0, AT_);
8310 MFloat distFac[2] = {18.0, 9.0};
8311 for(MInt i = 0; i < 2; i++) {
8312 distFac[i] = Context::getSolverProperty<MFloat>("mbBandWidth
", m_solverId, AT_, &distFac[i], i);
8315 m_outerBandWidth[grid().maxRefinementLevel() - 1] = distFac[0] * c_cellLengthAtLevel(grid().maxRefinementLevel());
8316 m_bandWidth[grid().maxRefinementLevel() - 1] = distFac[0];
8317 for(MInt i = grid().maxRefinementLevel() - 2; i >= 0; i--) {
8318 m_outerBandWidth[i] = m_outerBandWidth[i + 1] + (distFac[1] * c_cellLengthAtLevel(i + 1));
8319 m_bandWidth[i] = (m_bandWidth[i + 1] / 2) + 1 + distFac[1];
8321 for(MInt i = 0; i < grid().maxRefinementLevel(); i++) {
8322 m_innerBandWidth[i] = -m_outerBandWidth[i];
8323 m_log << "bandwidth level
" << i << ":
" << m_innerBandWidth[i] << " " << m_outerBandWidth[i] << endl;
8337 m_refineDiagonals = true;
8338 m_refineDiagonals = Context::getSolverProperty<MBool>("refineDiagonals
", m_solverId, AT_, &m_refineDiagonals);
8340 if(noNeighborDomains() > 0 && grid().isActive()) {
8341 grid().updateLeafCellExchange();
8349void LbSolver<nDim>::balancePre() {
8350 // Update the grid proxy for this solver
8353 // Just reset parallelization information if solver is not active and cleanup containers
8355 updateDomainInfo(-1, -1, MPI_COMM_NULL, AT_);
8359 // Set new domain info for solver
8360 updateDomainInfo(grid().domainId(), grid().noDomains(), grid().mpiComm(), AT_);
8362 // Clear cell collector
8365 // Resize cell collector to internal cells
8366 m_cells.append(grid().noInternalCells());
8368 // Now we are for the cell collector to be filled
8369 m_setCellDataFinished = false;
8376void LbSolver<nDim>::balancePost() {
8377 // Cell collector data has been set
8378 m_setCellDataFinished = true;
8380 // Nothing to do if solver is not active
8385 // append the halo-cells
8386 m_cells.append(c_noCells() - m_cells.size());
8388 updateCellCollectorFromGrid();
8390 grid().findEqualLevelNeighborsParDiagonal(false);
8394 // Restart boundary condition object to update cell lists
8396 resetActiveCellList();
8401 prepareCommunicationReduced();
8403 prepareCommunicationNormal();
8406 // Do an exchange to have the right values in halo cells
8407 if(noNeighborDomains() > 0) {
8409 exchangeOldDistributions();
8412 // If grid adaptation is used make sure that a grid file is written with the next restart file
8413 // since the cells will be sorted during balacing, while the old grid after adaptation will
8414 // contain the unsorted grid cells
8416 m_adaptationSinceLastRestart = true;
8417 m_adaptationSinceLastSolution = true;
8427MInt LbSolver<nDim>::noLoadTypes() const {
8428 const MInt noLevels = maxLevel() - minLevel() + 1;
8429 const MInt noBndryTypes = 1;
8431 return noLevels + noBndryTypes;
8440void LbSolver<nDim>::getDefaultWeights(MFloat* weights, std::vector<MString>& names) const {
8443 std::cerr << "WARNING: Using powers of 0.5 as weights. Could lead to large number of cells
for cases with large
"
8444 "refinement level range.
"
8449 const MInt noLevels = maxLevel() - minLevel() + 1;
8450 for(MInt level = 0; level < noLevels; level++) {
8451 const MInt lvlDiff = noLevels - level;
8452 weights[count] = FFPOW2(lvlDiff);
8453 std::stringstream ss;
8454 ss << "lb_leaf_cell_level_
" << minLevel() + level;
8455 names[count] = ss.str();
8460 weights[count] = 1.0;
8461 names[count] = "lb_bndry_cell
";
8464/* \brief Return the load of a single cell
8465 * \author Miro Gondrum
8467 * \param[in] cellId Requested grid cell id
8468 * \param[in] weights Computational weights for different simulation components
8472MFloat LbSolver<nDim>::getCellLoad(const MInt gridCellId, const MFloat* const weights) const {
8474 ASSERT(isActive(), "solver is not active
");
8476 // Convert to solver cell id and check
8477 const MInt cellId = grid().tree().grid2solver(gridCellId);
8482 if(cellId < 0 || cellId >= grid().noInternalCells()) {
8483 TERMM(1, "The given cell
id is
invalid.
");
8487 MFloat cellLoad = 0.0; // Default cell load
8488 if(a_isActive(cellId)) {
8489 const MInt relLevel = c_level(cellId) - minLevel();
8490 cellLoad += weights[relLevel];
8493 const MInt noLevels = maxLevel() - minLevel() + 1;
8494 if(a_isBndryCell(cellId)) {
8495 cellLoad += weights[noLevels];
8501// Get total number of loadQuantities for this domain...
8503void LbSolver<nDim>::getLoadQuantities(MInt* const loadQuantities) const {
8507 std::fill_n(&loadQuantities[0], noLoadTypes(), 0);
8509 // Nothing to do if solver is not active
8514 const MInt noLevels = maxLevel() - minLevel() + 1;
8516 for(MInt level = 0; level < noLevels; level++) {
8517 const MInt noActiveCellsOfLevel = m_activeCellListLvlOffset[level + 1] - m_activeCellListLvlOffset[level];
8518 loadQuantities[level] = noActiveCellsOfLevel;
8521 MInt noBoundaryCells = 0.0;
8522 for(MInt i = 0; i < a_noCells(); i++) {
8523 noBoundaryCells += a_isBndryCell(i);
8525 loadQuantities[noLevels] = noBoundaryCells;
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
MBool mDeallocate(T *&a)
deallocates the memory previously allocated for element 'a'
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
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(...).
Storage of the Position of the Conservative Variables (RHO, RHO_VV, RHO_E) in the value vectors of th...
Storage of the Position of the Primitive Variables (u, v, w, T, p) in the value vectors of the solver...
This class is a ScratchSpace.
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
@ QUADRATIC_INTERPOLATION
@ LB_TURBULENCE_ISOTROPIC_INIT
@ MAIA_LATTICE_BGK_TRANSPORT
@ MAIA_LATTICE_BGKI_EULER_2D
@ MAIA_LATTICE_BGKI_SMAGORINSKY
@ MAIA_LATTICE_BGK_TOTALENERGY_TRANSPORT
@ MAIA_LATTICE_BGK_TOTALENERGY
@ MAIA_LATTICE_BGK_INNERENERGY
@ MAIA_LATTICE_BGK_THERMAL_TRANSPORT
@ MAIA_LATTICE_RBGK_SMAGORINSKY
@ MAIA_LATTICE_CLB_SMAGORINSKY
@ MAIA_LATTICE_RBGK_DYNAMIC_SMAGO
@ MAIA_LATTICE_BGKI_SMAGORINSKY2
@ MAIA_LATTICE_BGK_INNERENERGY_TRANSPORT
@ MAIA_LATTICE_BGKI_DYNAMIC_SMAGO
@ MAIA_LATTICE_MRT_SMAGORINSKY
@ MAIA_LATTICE_BGK_THERMAL
@ MAIA_LATTICE_BGKI_SMAGO_WALL
@ MAIA_LATTICE_BGK_GUO_FORCING
void mTerm(const MInt errorCode, const MString &location, const MString &message)
constexpr Real POW3(const Real x)
constexpr Real POW2(const Real x)
double doubleSwap(double f)
constexpr MLong IPOW2(MInt x)
constexpr MFloat FPOW2(MInt x)
std::basic_string< char > MString
int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Status *status, const MString &name, const MString &varname)
same as MPI_Recv
int MPI_Comm_create(MPI_Comm comm, MPI_Group group, MPI_Comm *newcomm, const MString &name, const MString &varname)
same as MPI_Comm_create, but updates the number of MPI communicators
int MPI_Issend(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_Issend
int MPI_Group_incl(MPI_Group group, int n, const int ranks[], MPI_Group *newgroup, const MString &name)
same as MPI_Group_incl
int MPI_Comm_group(MPI_Comm comm, MPI_Group *group, const MString &name, const MString &varname)
same as MPI_Comm_group
int MPI_Wait(MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Wait
int MPI_Gatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, const int recvcounts[], const int displs[], MPI_Datatype recvtype, int root, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Gatherv
int MPI_Waitall(int count, MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Waitall
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_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_Gather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, int root, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Gather
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
IdType index(const FloatType *const x, const IdType level)
Return Hilbert index for given location and level in 2D or 3D.
MFloat determinant(std::array< T, N > &m)
MFloat distance(const MFloat *a, const MFloat *b)
MFloat frobeniusMatrixNormSquared(MFloatScratchSpace &m, MInt dim1, MInt dim2)
void addMatrices(MFloatScratchSpace &m1, MFloatScratchSpace &m2, MFloatScratchSpace &result, MInt dim1, MInt dim2)
void multiplyMatricesSq(MFloatScratchSpace &m1, MFloatScratchSpace &m2, MFloatScratchSpace &result, MInt dim)
void calcEigenValues(MFloat A[3][3], MFloat w[3])
void computeEnergySpectrum(MFloatScratchSpace &velocity, MIntScratchSpace &indices, const ptrdiff_t howmany, MInt locSize, MInt nx, MInt ny, MInt nz, MFloat viscosity, const MPI_Comm comm)
Compute energy spectrum on unity cube.
void exchangeData(const MInt noNghbrDomains, const MInt *const nghbrDomains, const MInt *const noHaloCells, const MInt **const, const MInt *const noWindowCells, const MInt **const windowCells, const MPI_Comm comm, const U *const data, U *const haloBuffer, const MInt noDat=1)
Generic exchange of data.
Namespace for auxiliary functions/classes.
void parallelFor(MInt begin, MInt end, UnaryFunction &&f)
Wrapper function for parallel for loop.
PARALLELIO_DEFAULT_BACKEND ParallelIo
MFloat dist(const Point< DIM > &p, const Point< DIM > &q)
LB lattice descriptor for arrays depending on D.
static constexpr MInt idFld(MInt i, MInt j)
static constexpr MInt oppositeDist(MInt i)
static constexpr MInt intNghbrArray(MInt i, MInt j)
static constexpr MFloat linearInterpolationCoefficients(MInt i, MInt j)
static constexpr MInt dirFld(MInt i, MInt j, MInt k)