27 m_grid->computeCellCenterCoordinates();
45 m_ransSolver = make_unique<FvStructuredSolver3DRans>(
this);
50 m_grid->computeCellCenterCoordinates();
63 if(
m_cells->
cellJac[cellId] < 0.0)
mTerm(1,
"Negative Jacobian found! Check first if code can cope this!");
122 cout <<
"///////////////////////////////////////////////////////////////////" << endl
123 <<
"Starting zonal bc creation for " << noZonalBCMaps <<
" zonal bc maps" << endl;
126 for(
MInt id = 0;
id < noZonalBCMaps;
id++) {
128 cout <<
"///////////////////////////////////////////////////////////////////" << endl;
129 cout <<
"//////////// ZONAL BC " <<
id <<
" ///////////////////////////////////////////" << endl;
130 cout <<
"///////////////////////////////////////////////////////////////////" << endl;
137 cout <<
"Zonal BC" <<
id <<
": Initialization ..." << endl;
142 MInt hasRcvDomain = 0;
162 m_zonalBC[
id]->noCellsLocalBC += (end[0] + stgIP - start[0]) * (end[1] - start[1]) * (end[2] - start[2]);
175 cout <<
"Zonal BC" <<
id <<
": Initialization ... Finished!" << endl;
176 cout <<
"Zonal BC" <<
id <<
": Collecting local coordinates ... " << endl;
191 localMapCellIds[bcCellId] = bcCellId;
193 localMapCellIds[bcCellId] = cellId;
196 for(
MInt dim = 0; dim <
nDim; ++dim) {
200 localReceiverIds[bcCellId] =
domainId();
208 cout <<
"Zonal BC" <<
id <<
": Collecting local coordinates ... Finished!" << endl;
209 cout <<
"Zonal BC" <<
id <<
": Exchanging information about local coordinates ..." << endl;
213 m_StructuredComm, AT_,
"m_zonalBC[id]->noCellsLocalBC",
"m_zonalBC[id]->noCellsGlobalBC");
215 "m_zonalBC[id]->noCellsLocalBC",
"noZonalCells[0]");
218 cout <<
"Zonal BC" <<
id <<
": Exchanging information about local coordinates ... Finished!" << endl;
227 noCells += noZonalCells[i];
229 cout <<
"Zonal BC" <<
id <<
": Gathering coordinates for " << noCells <<
" cells on all partitions..." << endl;
237 mAlloc(
m_zonalBC[
id]->globalReceiverIds, std::max(1,
m_zonalBC[
id]->noCellsGlobalBC),
"globalReceiverIds", -1, AT_);
239 noCellsBCOffsets[0] = 0;
241 noCellsBCOffsets[i] = noCellsBCOffsets[i - 1] + noZonalCells[i - 1];
246 for(
MInt dim = 0; dim <
nDim; ++dim) {
248 &
m_zonalBC[
id]->coordinatesGlobalBC[dim][0], &noZonalCells[0], &noCellsBCOffsets[0], MPI_DOUBLE,
249 m_StructuredComm, AT_,
"localCoordinatesBC(0,0)",
"coordinatesGlobalBC[0][0]");
252 &noZonalCells[0], &noCellsBCOffsets[0], MPI_INT,
m_StructuredComm, AT_,
"localReceiverIds[0]",
253 "m_zonalBC[id]->globalReceiverIds[0]");
255 &noZonalCells[0], &noCellsBCOffsets[0], MPI_INT,
m_StructuredComm, AT_,
"localMapCellIds[0]",
256 "m_zonalBC[id]->globalLocalMapCellIds[0]");
261 noCells += noZonalCells[i];
263 cout <<
"Zonal BC" <<
id <<
": Gathering coordinates for " << noCells <<
" cells on all partitions... Finished!"
265 cout <<
"Zonal BC" <<
id <<
": Searching for suitable interpolation partners..." << endl;
271 MBool hasInterpolationPartnerDomain =
true;
274 hasInterpolationPartnerDomain =
false;
282 hasPartnerLocalBC.
begin(),
283 hasInterpolationPartnerDomain);
286 cout <<
"Zonal BC" <<
id <<
": Searching for suitable interpolation partners... Finished!" << endl;
287 cout <<
"Zonal BC" <<
id <<
": Creating sending and receiving information..." << endl;
297 "hasRcvDomain",
"m_zonalBC[id]->noGlobalRcvDomains");
304 if(noZonalCells[j] > 0) {
311 MInt hasSndDomain =
false;
312 for(
MInt cellIdBC = 0; cellIdBC <
m_zonalBC[
id]->noCellsGlobalBC; cellIdBC++) {
313 if(hasPartnerLocalBC[cellIdBC]) {
321 "hasSndDomain",
"m_zonalBC[id]->noGlobalSndDomains");
326 "hasSndDomainInfo[0]");
333 if(hasSndDomainInfo[j] > 0) {
343 for(
MInt cellIdBC = 0; cellIdBC <
m_zonalBC[
id]->noCellsGlobalBC; cellIdBC++) {
344 if(hasPartnerLocalBC[cellIdBC]) {
358 for(
MInt cellIdBC = 0; cellIdBC <
m_zonalBC[
id]->noCellsGlobalBC; cellIdBC++) {
359 if(hasPartnerLocalBC[cellIdBC]) {
371 if(
m_zonalBC[
id]->localCommReceiverIds[j] == d) {
389 if(
m_zonalBC[
id]->localCommReceiverIds[j] == d) {
402 for(
MInt cellIdBC = 0; cellIdBC <
m_zonalBC[
id]->noCellsGlobalBC; cellIdBC++) {
403 if(hasPartnerLocalBC[cellIdBC]) {
405 m_zonalBC[
id]->localBufferMapCellIds[pos] = mapCellId;
406 m_zonalBC[
id]->localBufferIndexCellIds[pos] = cellIdBC;
414 if(localRcvId[i] ==
m_zonalBC[
id]->localCommReceiverIds[j]) {
415 localBufferSndSize[i]++;
423 for(
MInt noRcv = 0; noRcv <
m_zonalBC[
id]->noGlobalRcvDomains; noRcv++) {
424 MInt bufferRcvSize = 0;
425 for(
MInt cellIdBC = 0; cellIdBC <
m_zonalBC[
id]->noCellsGlobalBC; cellIdBC++) {
426 if(
m_zonalBC[
id]->globalRcvDomainIds[noRcv] ==
m_zonalBC[
id]->globalReceiverIds[cellIdBC]) {
427 if(hasPartnerLocalBC[cellIdBC]) {
434 sndBufferRcvSize[
m_zonalBC[
id]->globalRcvDomainIds[noRcv]] = bufferRcvSize;
439 "sndBufferRcvSize[0]",
"rcvBufferRcvSize[0]");
447 if(rcvBufferRcvSize[i] > 0) {
460 if(rcvBufferRcvSize[i] > 0) {
470 if(rcvBufferRcvSize[i] > 0) {
471 localBufferRcvSize[pos] = rcvBufferRcvSize[i];
488 unique_ptr<StructuredZonalComm> sndCommPtr =
489 make_unique<StructuredZonalComm>(localBufferSndSize[i], localRcvId[i],
m_zonalBC[
id]->noZonalVariables);
490 m_zonalBC[
id]->sndComm.push_back(std::move(sndCommPtr));
495 unique_ptr<StructuredZonalComm> rcvCommPtr =
496 make_unique<StructuredZonalComm>(localBufferRcvSize[i], localSndId[i],
m_zonalBC[
id]->noZonalVariables);
497 m_zonalBC[
id]->rcvComm.push_back(std::move(rcvCommPtr));
506 for(
MInt noSnd = 0; noSnd <
m_zonalBC[
id]->noSndNghbrDomains; noSnd++) {
509 for(
MInt cellId = 0; cellId <
m_zonalBC[
id]->noBufferSndSize; cellId++) {
510 if(
m_zonalBC[
id]->localCommReceiverIds[cellId] ==
m_zonalBC[
id]->sndComm[noSnd]->localId) {
518 &
m_zonalBC[
id]->sndRequest[noSnd], AT_,
"m_zonalBC[id]->sndComm[noSnd]->intBbuffer");
520 cout <<
"rank " <<
domainId() <<
" zonal sending to " <<
m_zonalBC[
id]->sndComm[noSnd]->localId
521 <<
" throws error " << endl;
524 for(
MInt noRcv = 0; noRcv <
m_zonalBC[
id]->noRcvNghbrDomains; noRcv++) {
528 &
m_zonalBC[
id]->rcvRequest[noRcv], AT_,
"(void*)&m_zonalBC[id]->rcvComm[noRcv]->mapCellId");
530 cout <<
"rank " <<
domainId() <<
" zonal sending to " <<
m_zonalBC[
id]->sndComm[noRcv]->localId
531 <<
" throws error " << endl;
540 cout <<
"Zonal BC" <<
id <<
": Creating sending and receiving information... Finished!" << endl;
541 cout <<
"Zonal BC" <<
id <<
": Completed!" << endl;
557 for(
MInt noSnd = 0; noSnd < zBC->noSndNghbrDomains; noSnd++) {
558 const MInt noCells = zBC->sndComm[noSnd]->bufferSize;
560 for(
MInt cellId = 0; cellId < zBC->noBufferSndSize; cellId++) {
561 if(zBC->localCommReceiverIds[cellId] == zBC->sndComm[noSnd]->localId) {
562 const MInt cellIdBC = zBC->localBufferIndexCellIds[cellId];
564 zBC->sndComm[noSnd]->buffer[pos + 0 * noCells] =
566 zBC->sndComm[noSnd]->buffer[pos + 1 * noCells] =
568 zBC->sndComm[noSnd]->buffer[pos + 2 * noCells] =
570 zBC->sndComm[noSnd]->buffer[pos + 3 * noCells] =
572 zBC->sndComm[noSnd]->buffer[pos + 4 * noCells] =
574 zBC->sndComm[noSnd]->buffer[pos + 5 * noCells] =
575 zBC->interpolation->interpolateVariableZonal(
m_cells->
fq[
FQ->NU_T], cellIdBC);
577 zBC->sndComm[noSnd]->buffer[pos + 0 * noCells] =
578 zBC->interpolation->interpolateVariableZonal(
m_cells->
fq[
FQ->AVG_U], cellIdBC);
579 zBC->sndComm[noSnd]->buffer[pos + 1 * noCells] =
580 zBC->interpolation->interpolateVariableZonal(
m_cells->
fq[
FQ->AVG_V], cellIdBC);
581 zBC->sndComm[noSnd]->buffer[pos + 2 * noCells] =
582 zBC->interpolation->interpolateVariableZonal(
m_cells->
fq[
FQ->AVG_W], cellIdBC);
583 zBC->sndComm[noSnd]->buffer[pos + 3 * noCells] =
584 zBC->interpolation->interpolateVariableZonal(
m_cells->
fq[
FQ->AVG_RHO], cellIdBC);
585 zBC->sndComm[noSnd]->buffer[pos + 4 * noCells] =
586 zBC->interpolation->interpolateVariableZonal(
m_cells->
fq[
FQ->AVG_P], cellIdBC);
587 zBC->sndComm[noSnd]->buffer[pos + 5 * noCells] =
588 zBC->interpolation->interpolateVariableZonal(
m_cells->
fq[
FQ->NU_T], cellIdBC);
600 for(
MInt i = 0; i < zBC->noGlobalSndDomains; i++) {
601 if(zBC->globalSndDomainIds[i] ==
domainId()) {
602 for(
MInt noSnd = 0; noSnd < zBC->noSndNghbrDomains; noSnd++) {
603 MInt err =
MPI_Isend(&zBC->sndComm[noSnd]->buffer[0], zBC->sndComm[noSnd]->bufferSize * zBC->noZonalVariables,
604 MPI_DOUBLE, zBC->sndComm[noSnd]->localId, 0,
m_StructuredComm, &zBC->sndRequest[noSnd],
605 AT_,
"zBC->buffer[noSnd][0]");
606 if(err) cout <<
"rank " <<
domainId() <<
" zonal sending throws error " << endl;
615 for(
MInt i = 0; i < zBC->noGlobalRcvDomains; i++) {
617 if(zBC->globalRcvDomainIds[i] ==
domainId()) {
618 for(
MInt noRcv = 0; noRcv < zBC->noRcvNghbrDomains; noRcv++) {
619 MInt err =
MPI_Irecv(&zBC->rcvComm[noRcv]->buffer[0], zBC->rcvComm[noRcv]->bufferSize * zBC->noZonalVariables,
620 MPI_DOUBLE, zBC->rcvComm[noRcv]->localId, 0,
m_StructuredComm, &zBC->rcvRequest[noRcv],
621 AT_,
"&zBC->bufferRcvZonal[noRcv][0]");
622 if(err) cout <<
"rank " <<
domainId() <<
" zonal receiving throws error " << endl;
626 MPI_Waitall(zBC->noSndNghbrDomains, zBC->sndRequest, zBC->sndStatus, AT_);
627 MPI_Waitall(zBC->noRcvNghbrDomains, zBC->rcvRequest, zBC->rcvStatus, AT_);
634 for(
MInt noRcv = 0; noRcv < zBC->noRcvNghbrDomains; noRcv++) {
635 const MInt noCells = zBC->rcvComm[noRcv]->bufferSize;
636 for(
MInt i = 0; i < noCells; i++) {
637 const MInt localMapCellIds = zBC->rcvComm[noRcv]->mapCellId[i];
639 m_cells->
stg_fq[
PV->U][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[0 * noCells + i];
640 m_cells->
stg_fq[
PV->V][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[1 * noCells + i];
641 m_cells->
stg_fq[
PV->W][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[2 * noCells + i];
642 m_cells->
stg_fq[
PV->RHO][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[3 * noCells + i];
643 m_cells->
stg_fq[
PV->P][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[4 * noCells + i];
644 m_cells->
stg_fq[
FQ->NU_T][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[5 * noCells + i];
646 m_cells->
fq[
FQ->AVG_U][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[0 * noCells + i];
647 m_cells->
fq[
FQ->AVG_V][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[1 * noCells + i];
648 m_cells->
fq[
FQ->AVG_W][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[2 * noCells + i];
649 m_cells->
fq[
FQ->AVG_RHO][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[3 * noCells + i];
650 m_cells->
fq[
FQ->AVG_P][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[4 * noCells + i];
651 m_cells->
fq[
FQ->NU_T][localMapCellIds] = zBC->rcvComm[noRcv]->buffer[5 * noCells + i];
689 m_log <<
"Using VENKATAKRISHNAN MOD limiter with VENK factor of " <<
m_venkFactor <<
" !" << endl;
696 m_log <<
"Using VENKATAKRISHNAN limiter!" << endl;
703 m_log <<
"Using BARTH JESPERSON limiter!" << endl;
709 m_log <<
"Using MINMOD limiter!" << endl;
714 m_log <<
"Using VAN ALBADA limiter!" << endl;
719 stringstream errorMessage;
720 errorMessage <<
"Limiter function " <<
m_limiterMethod <<
" not implemented!" << endl;
721 mTerm(1, AT_, errorMessage.str());
725 m_log <<
"Using unlimited MUSCL! (standard Formulation)" << endl;
727 m_log <<
"Using standard AUSM central" << endl;
729 switch(
CV->noVariables) {
743 stringstream errorMessage;
744 errorMessage <<
"Number of Variables " <<
CV->noVariables <<
" not implemented in template AUSM!" << endl;
749 m_log <<
"Using AUSM PTHRC" << endl;
751 switch(
CV->noVariables) {
765 stringstream errorMessage;
766 errorMessage <<
"Number of Variables " <<
CV->noVariables <<
" not implemented in template AUSM!" << endl;
771 m_log <<
"Using AUSMDV" << endl;
773 switch(
CV->noVariables) {
787 stringstream errorMessage;
788 errorMessage <<
"Number of Variables " <<
CV->noVariables <<
" not implemented in template AUSM!" << endl;
794 m_log <<
"Using unlimited MUSCL (streched Grids)";
797 m_log <<
"Using standard AUSM central" << endl;
799 switch(
CV->noVariables) {
813 stringstream errorMessage;
814 errorMessage <<
"Number of Variables " <<
CV->noVariables <<
" not implemented in template AUSM!" << endl;
819 m_log <<
"Using AUSM PTHRC" << endl;
820 switch(
CV->noVariables) {
834 stringstream errorMessage;
835 errorMessage <<
"Number of Variables " <<
CV->noVariables <<
" not implemented in template AUSM!" << endl;
840 m_log <<
"Using AUSMDV" << endl;
841 switch(
CV->noVariables) {
855 stringstream errorMessage;
856 errorMessage <<
"Number of Variables " <<
CV->noVariables <<
" not implemented in template AUSM!" << endl;
870 MInt P1 = -1, P2 = -1, P3 = -1, P4 = -1, P5 = -1, P6 = -1, P7 = -1, P8 = -1;
871 MFloat f1x = F0, f1y = F0, f1z = F0, f2x = F0, f2y = F0, f2z = F0;
887 * (
m_grid->m_coordinates[0][P1] +
m_grid->m_coordinates[0][P6] +
m_grid->m_coordinates[0][P7]
888 +
m_grid->m_coordinates[0][P8]);
890 * (
m_grid->m_coordinates[1][P1] +
m_grid->m_coordinates[1][P6] +
m_grid->m_coordinates[1][P7]
891 +
m_grid->m_coordinates[1][P8]);
893 * (
m_grid->m_coordinates[2][P1] +
m_grid->m_coordinates[2][P6] +
m_grid->m_coordinates[2][P7]
894 +
m_grid->m_coordinates[2][P8]);
897 * (
m_grid->m_coordinates[0][P2] +
m_grid->m_coordinates[0][P3] +
m_grid->m_coordinates[0][P4]
898 +
m_grid->m_coordinates[0][P5]);
900 * (
m_grid->m_coordinates[1][P2] +
m_grid->m_coordinates[1][P3] +
m_grid->m_coordinates[1][P4]
901 +
m_grid->m_coordinates[1][P5]);
903 * (
m_grid->m_coordinates[2][P2] +
m_grid->m_coordinates[2][P3] +
m_grid->m_coordinates[2][P4]
904 +
m_grid->m_coordinates[2][P5]);
909 * (
m_grid->m_coordinates[0][P1] +
m_grid->m_coordinates[0][P2] +
m_grid->m_coordinates[0][P4]
910 +
m_grid->m_coordinates[0][P8]);
912 * (
m_grid->m_coordinates[1][P1] +
m_grid->m_coordinates[1][P2] +
m_grid->m_coordinates[1][P4]
913 +
m_grid->m_coordinates[1][P8]);
915 * (
m_grid->m_coordinates[2][P1] +
m_grid->m_coordinates[2][P2] +
m_grid->m_coordinates[2][P4]
916 +
m_grid->m_coordinates[2][P8]);
919 * (
m_grid->m_coordinates[0][P3] +
m_grid->m_coordinates[0][P5] +
m_grid->m_coordinates[0][P6]
920 +
m_grid->m_coordinates[0][P7]);
922 * (
m_grid->m_coordinates[1][P3] +
m_grid->m_coordinates[1][P5] +
m_grid->m_coordinates[1][P6]
923 +
m_grid->m_coordinates[1][P7]);
925 * (
m_grid->m_coordinates[2][P3] +
m_grid->m_coordinates[2][P5] +
m_grid->m_coordinates[2][P6]
926 +
m_grid->m_coordinates[2][P7]);
931 * (
m_grid->m_coordinates[0][P1] +
m_grid->m_coordinates[0][P2] +
m_grid->m_coordinates[0][P3]
932 +
m_grid->m_coordinates[0][P6]);
934 * (
m_grid->m_coordinates[1][P1] +
m_grid->m_coordinates[1][P2] +
m_grid->m_coordinates[1][P3]
935 +
m_grid->m_coordinates[1][P6]);
937 * (
m_grid->m_coordinates[2][P1] +
m_grid->m_coordinates[2][P2] +
m_grid->m_coordinates[2][P3]
938 +
m_grid->m_coordinates[2][P6]);
941 * (
m_grid->m_coordinates[0][P4] +
m_grid->m_coordinates[0][P5] +
m_grid->m_coordinates[0][P7]
942 +
m_grid->m_coordinates[0][P8]);
944 * (
m_grid->m_coordinates[1][P4] +
m_grid->m_coordinates[1][P5] +
m_grid->m_coordinates[1][P7]
945 +
m_grid->m_coordinates[1][P8]);
947 * (
m_grid->m_coordinates[2][P5] +
m_grid->m_coordinates[2][P5] +
m_grid->m_coordinates[2][P7]
948 +
m_grid->m_coordinates[2][P8]);
965 const MFloat lowerY = F1B4
966 * (
m_grid->m_coordinates[1][P1] +
m_grid->m_coordinates[1][P2] +
m_grid->m_coordinates[1][P4]
967 +
m_grid->m_coordinates[1][P8]);
968 const MFloat upperY = F1B4
969 * (
m_grid->m_coordinates[1][P3] +
m_grid->m_coordinates[1][P5] +
m_grid->m_coordinates[1][P6]
970 +
m_grid->m_coordinates[1][P7]);
972 return upperY - lowerY;
1069 MBool donorConservative =
false;
1071 donorConservative = Context::getSolverProperty<MBool>(
"donorConservative",
false, AT_);
1074 if(donorConservative) {
1075 for(
MInt var = 0; var <
CV->noVariables; var++) {
1080 for(
MInt var = 0; var <
PV->noVariables; var++) {
1118 PV->TInfinity = 1.0 / (1.0 + F1B2 * gammaMinusOne *
POW2(
m_Ma));
1119 UT =
m_Ma * sqrt(
PV->TInfinity);
1123 PV->VVInfinity[0] =
PV->UInfinity;
1124 PV->VVInfinity[1] =
PV->VInfinity;
1125 PV->VVInfinity[2] =
PV->WInfinity;
1129 CV->rhoInfinity = pow(
PV->TInfinity, (1.0 / gammaMinusOne));
1130 CV->rhoUInfinity =
CV->rhoInfinity *
PV->UInfinity;
1131 CV->rhoVInfinity =
CV->rhoInfinity *
PV->VInfinity;
1132 CV->rhoWInfinity =
CV->rhoInfinity *
PV->WInfinity;
1133 CV->rhoVVInfinity[0] =
CV->rhoUInfinity;
1134 CV->rhoVVInfinity[1] =
CV->rhoVInfinity;
1135 CV->rhoVVInfinity[2] =
CV->rhoWInfinity;
1136 CV->rhoEInfinity =
PV->PInfinity / gammaMinusOne +
CV->rhoInfinity * (F1B2 *
POW2(UT));
1139 m_Re0 =
m_Re * SUTHERLANDLAW(
PV->TInfinity) / (
CV->rhoInfinity *
m_Ma * sqrt(
PV->TInfinity));
1155 const MFloat lamVisc = SUTHERLANDLAW(
PV->TInfinity);
1157 CV->ransInfinity[0] = chi * (lamVisc);
1158 PV->ransInfinity[0] = chi * (lamVisc /
CV->rhoInfinity);
1161 m_log <<
"=================================================" << endl;
1162 m_log <<
" INITIAL CONDITION SUMMARY" << endl;
1163 m_log <<
"=================================================" << endl;
1167 m_log <<
"TInfinity = " <<
PV->TInfinity << endl;
1168 m_log <<
"UInfinity = " <<
PV->UInfinity << endl;
1169 m_log <<
"VInfinity = " <<
PV->VInfinity << endl;
1170 m_log <<
"WInfinity = " <<
PV->WInfinity << endl;
1171 m_log <<
"PInfinity = " <<
PV->PInfinity << endl;
1172 m_log <<
"rhoInfinity = " <<
CV->rhoInfinity << endl;
1173 m_log <<
"rhoEInfinity = " <<
CV->rhoEInfinity << endl;
1177 cout <<
"////////////////////////////////////////////////" << endl;
1178 cout <<
"////////// Initial Condition summary ///////////" << endl;
1179 cout <<
"////////////////////////////////////////////////" << endl;
1180 cout <<
"Re = " <<
m_Re << endl;
1181 cout <<
"Re0 = " <<
m_Re0 << endl;
1182 cout <<
"Ma = " <<
m_Ma << endl;
1183 cout <<
"TInfinity = " <<
PV->TInfinity << endl;
1184 cout <<
"UInfinity = " <<
PV->UInfinity << endl;
1185 cout <<
"VInfinity = " <<
PV->VInfinity << endl;
1186 cout <<
"WInfinity = " <<
PV->WInfinity << endl;
1188 cout <<
"PInfinity = " <<
PV->PInfinity << endl;
1189 cout <<
"rhoInfinity = " <<
CV->rhoInfinity << endl;
1190 cout <<
"rhoEInfinity = " <<
CV->rhoEInfinity << endl;
1191 cout <<
"referenceTime = " <<
m_timeRef << endl;
1192 cout <<
" zonal = " <<
m_zonal << endl;
1232 if(radius <= 0.05) {
1234 MFloat pressureSignal = sin(radius / 0.05 * PI) * pAmp +
PV->PInfinity;
1255 MFloat pressureSignal = sin(xCoordinate / 0.1 * PI) * pAmp +
PV->PInfinity;
1271 cout <<
"I.C. stagnating flow field was applied! " << endl;
1288 -(F3 / F2) *
PV->UInfinity * (
POW2(
y - y_max / F2) -
POW2(y_max / F2)) /
POW2(y_max / F2);
1292 PV->PInfinity - F3 * (x + 15.0) * SUTHERLANDLAW(
PV->TInfinity) *
PV->UInfinity *
POW2(F2 / y_max) /
m_Re0;
1331 + F1B16 * (
POW2(A) *
CV->rhoInfinity) * (cos(2.0 * x) + cos(2.0 *
y)) * (2.0 + cos(2.0 * z));
1397 m_log <<
"deltaP: " << deltaP <<
" cfTheoretisch: " << cfTheo <<
" cdTheo: " << cdTheo << endl;
1415 yplus = prefactor *
y;
1423 }
else if(yplus <= 30 && yplus > 5.0) {
1425 }
else if(yplus > 30) {
1445 const MInt totalBlockCells =
1446 (
m_grid->getMyBlockNoCells(0) *
m_grid->getMyBlockNoCells(1) *
m_grid->getMyBlockNoCells(2));
1448 if(totalBlockCells <= 12500000) {
1449 fftw_complex *uPhysField, *vPhysField, *wPhysField;
1452 MInt lx =
m_grid->getMyBlockNoCells(2), ly =
m_grid->getMyBlockNoCells(1), lz =
m_grid->getMyBlockNoCells(0);
1455 uPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
1456 vPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
1457 wPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
1467 MInt m_noPeakModes = 100;
1468 initFFTW(uPhysField, vPhysField, wPhysField, lx, ly, lz, m_noPeakModes);
1471 for(
MInt id = 0;
id < lz * ly * lz;
id++) {
1472 sendRecvBufferU[
id] = uPhysField[
id][0];
1473 sendRecvBufferV[
id] = vPhysField[
id][0];
1474 sendRecvBufferW[
id] = wPhysField[
id][0];
1482 for(
MInt id = 0;
id < lz * ly * lz;
id++) {
1483 uPhysField[
id][0] = sendRecvBufferU[
id];
1484 vPhysField[
id][0] = sendRecvBufferV[
id];
1485 wPhysField[
id][0] = sendRecvBufferW[
id];
1490 MInt m_noPeakModes = 100;
1491 initFFTW(uPhysField, vPhysField, wPhysField, lx, ly, lz, m_noPeakModes);
1501 amp * 2.0 * (0.5 - (1.0 * rand() / (RAND_MAX + 1.0))) *
m_cells->
pvariables[
PV->U][cellId];
1503 amp * 2.0 * (0.5 - (1.0 * rand() / (RAND_MAX + 1.0))) *
m_cells->
pvariables[
PV->U][cellId];
1505 amp * 2.0 * (0.5 - (1.0 * rand() / (RAND_MAX + 1.0))) *
m_cells->
pvariables[
PV->U][cellId];
1530 m_log <<
"=========== Turb. Pipe Flow Inital Condition Summary =========== " << endl;
1531 m_log <<
"-->Turbulent pipe flow deltaP: " <<
m_deltaP << endl;
1532 m_log <<
"-->pipe friciton velocity: " << uTau << endl;
1535 m_log <<
"--> bulk velocity (u_max)" << bulkVel << endl;
1536 m_log <<
"=========== Turb. Pipe Flow Initial Condition Summary Finished =========== " << endl;
1571 const MInt totalBlockCells =
1572 (
m_grid->getMyBlockNoCells(0) *
m_grid->getMyBlockNoCells(1) *
m_grid->getMyBlockNoCells(2));
1574 if(totalBlockCells <= 12500000) {
1575 fftw_complex *uPhysField, *vPhysField, *wPhysField;
1578 MInt lx =
m_grid->getMyBlockNoCells(2), ly =
m_grid->getMyBlockNoCells(1), lz =
m_grid->getMyBlockNoCells(0);
1581 uPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
1582 vPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
1583 wPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
1593 MInt m_noPeakModes = 100;
1594 initFFTW(uPhysField, vPhysField, wPhysField, lx, ly, lz, m_noPeakModes);
1597 for(
MInt id = 0;
id < lz * ly * lz;
id++) {
1598 sendRecvBufferU[
id] = uPhysField[
id][0];
1599 sendRecvBufferV[
id] = vPhysField[
id][0];
1600 sendRecvBufferW[
id] = wPhysField[
id][0];
1608 for(
MInt id = 0;
id < lz * ly * lz;
id++) {
1609 uPhysField[
id][0] = sendRecvBufferU[
id];
1610 vPhysField[
id][0] = sendRecvBufferV[
id];
1611 wPhysField[
id][0] = sendRecvBufferW[
id];
1616 MInt m_noPeakModes = 100;
1617 initFFTW(uPhysField, vPhysField, wPhysField, lx, ly, lz, m_noPeakModes);
1627 amp * 2.0 * (0.5 - (1.0 * rand() / (RAND_MAX + 1.0))) *
m_cells->
pvariables[
PV->U][cellId];
1629 amp * 2.0 * (0.5 - (1.0 * rand() / (RAND_MAX + 1.0))) *
m_cells->
pvariables[
PV->U][cellId];
1631 amp * 2.0 * (0.5 - (1.0 * rand() / (RAND_MAX + 1.0))) *
m_cells->
pvariables[
PV->U][cellId];
1641 const MFloat epss = 1e-10;
1642 const MFloat reTheta = 1000.000;
1643 const MFloat theta = 1.0;
1644 const MFloat delta0 = 72.0 / 7.0 * theta;
1645 const MFloat kappa = 0.4;
1646 const MFloat C1 = 3.573244189003983;
1648 const MFloat cf = 0.024 / pow(reTheta, 0.25);
1649 const MFloat uTau = sqrt(cf / 2.0) *
PV->UInfinity;
1650 const MFloat nu = SUTHERLANDLAW(
PV->TInfinity);
1661 }
else if(yPlus <= 10.0) {
1663 }
else if(yPlus <= 30 && yPlus > 10.0) {
1665 }
else if(yPlus > 30.0 &&
y <= delta0) {
1668 * ((F1 / kappa) * log(max(yPlus, epss)) + C1 + 2 * (PI1 / kappa) * (3 *
y *
y - 2 *
y *
y *
y));
1688 const MFloat epss = 1e-10;
1689 const MFloat reTheta = 1000.0;
1690 const MFloat theta = 1.0;
1691 const MFloat delta0 = 72.0 / 7.0 * theta;
1693 const MFloat C1 = 3.573244189003983;
1695 const MFloat cf = 0.024 / pow(reTheta, 0.25);
1701 const MFloat mu = SUTHERLANDLAW(
PV->TInfinity);
1702 const MFloat utau = sqrt(cf / 2.0) *
m_Ma * sqrt(
PV->TInfinity);
1710 }
else if(yplus < 10) {
1715 * ((1. / K) * log(max(yplus, epss)) + C1 + 2 * PI1 / K * (3 * eta * eta - 2 * eta * eta * eta)),
1751 MFloat disturb = amp * exp(-(a1 + a2 + a3) /
POW2(r) / 2.0);
1762 MFloat alpha = log(2) / 9;
1803 amp * (r0 / r) * (r - r0) * exp(-1.0 * alpha * (
POW2(x) +
POW2(r - r0)));
1804 v = -1.0 * amp * (r0 / r) * x * exp(-1.0 * alpha * (
POW2(x) +
POW2(r - r0)));
1831 amp * (r0 / r) * (r - r0) * exp(-1.0 * alpha * (
POW2(x) +
POW2(r - r0)));
1832 v = -1.0 * amp * (r0 / r) * x * exp(-1.0 * alpha * (
POW2(x) +
POW2(r - r0)));
1870 MFloat u = F1B2 * (F1 - tanh(12.5 * (fabs(r / 0.5) - fabs(0.5 / r)))) *
PV->VVInfinity[0];
1939 amp * (r0 /
mMax(r, 0.00001)) * (r - r0) * exp(-1.0 * alpha * (
POW2(x) +
POW2(r - r0)));
1940 v = -1.0 * amp * (r0 /
mMax(r, 0.00001)) * x * exp(-1.0 * alpha * (
POW2(x) +
POW2(r - r0)));
1948 m_log <<
"falkner skan cooke initial condition (incompressible)" << endl;
1949 if(!
m_fsc)
mTerm(1,
"property fsc not set. Refer to the description of the property");
1953 fscf.open(
"pressure.dat", ios::trunc);
1955 fscf <<
"#x_maia x_fsc p" << endl;
1956 for(
MInt i = 0; i < 95; i++) {
1963 fscf.open(
"velocity_x0.dat", ios::trunc);
1965 fscf <<
"#y eta u v w" << endl;
1966 for(
MInt i = 0; i < 200; i++) {
1973 for(
MInt dim = 0; dim <
nDim; dim++)
1974 fscf <<
" " << vel[dim] /
PV->UInfinity;
1996 m_log <<
"Blasius initial condition (incompressible)" << endl;
1997 if(!
m_useBlasius)
mTerm(1,
"property Blasius not set. Refer to the description of the property");
2001 blasiusf.open(
"velocity_x0.dat", ios::trunc);
2003 blasiusf <<
"#y eta u v" << endl;
2007 MBool d0Set =
false;
2016 for(
MInt dim = 0; dim <
nDim; dim++)
2017 blasiusf <<
" " << vel[dim] /
PV->UInfinity;
2020 if(!d0Set && vel[0] >= 0.99 *
PV->UInfinity) {
2026 d1 += (1 - vel[0] /
PV->UInfinity) * 0.05;
2027 d2 += vel[0] /
PV->UInfinity * (1 - vel[0] /
PV->UInfinity) * 0.05;
2033 cout <<
"d0: " << d0 <<
" d1: " << d1 <<
" d2: " << d2 << endl;
2052 mTerm(1, AT_,
"No (correct) initial Condition is given!");
2072 MPI_SUM,
m_StructuredComm, AT_,
"m_intpPointsHasPartnerLocal",
"m_intpPointsHasPartnerGlobal");
2087 MPI_SUM,
m_StructuredComm, AT_,
"m_pointsToAsciiHasPartnerLocal",
"m_pointsToAsciiHasPartnerGlobal");
2090 cout <<
"domainId: " <<
domainId() <<
" pointId: " << pointId
2113 for(
MInt var = 0; var <
PV->noVariables; var++) {
2126 "m_pointsToAsciiHasPartnerGlobal");
2139 MPI_SUM,
m_StructuredComm, AT_,
"m_pointsToAsciiVars[0][0]",
"m_pointsToAsciiVars[0][0]");
2159 cout <<
"globalTimeStep: " <<
globalTimeStep <<
" writing point data out to ascii file" << endl;
2160 MString filename =
"./pointVars.dat";
2162 f_forces = fopen(filename.c_str(),
"a+");
2171 fprintf(f_forces,
"\n");
2203 MBool interpolationCorrection =
false;
2205 interpolationCorrection =
2206 Context::getSolverProperty<MBool>(
"interpolationCorrection",
m_solverId, AT_, &interpolationCorrection);
2209 cout <<
"Correcting interpolation" << endl;
2216 for(
MInt var = 0; var <
PV->noVariables; var++) {
2235 static MBool deviateAvailable =
false;
2236 static float storedDeviate;
2237 MFloat polar, rsquared, var1, var2;
2242 if(!deviateAvailable) {
2246 var1 = 2.0 * (
MFloat(rand()) /
MFloat(RAND_MAX)) - 1.0;
2247 var2 = 2.0 * (
MFloat(rand()) /
MFloat(RAND_MAX)) - 1.0;
2248 rsquared = var1 * var1 + var2 * var2;
2249 }
while(rsquared >= 1.0 ||
approx(rsquared, F0,
m_eps));
2252 polar = sqrt(-2.0 * log(rsquared) / rsquared);
2255 storedDeviate = var1 * polar;
2256 deviateAvailable =
true;
2259 return var2 * polar * sigma + mu;
2264 deviateAvailable =
false;
2265 return storedDeviate * sigma + mu;
2289 const MInt sword = 2;
2300 for(
MInt dim = 0; dim <
nDim; ++dim) {
2306 const MInt IP1 = I + IJK[dim];
2307 const MInt IM1 = I - IJK[dim];
2308 const MInt IP2 = I + 2 * IJK[dim];
2314 const MFloat DSP1 = sqrt(
POW2(x[IP2] - x[IP1]) +
POW2(
y[IP2] -
y[IP1]) +
POW2(z[IP2] - z[IP1]));
2319 for(
MInt var = 0; var <
PV->noVariables; ++var) {
2320 const MUint offset = var * noCells;
2321 const MFloat*
const RESTRICT vars = ALIGNED_F(cellVariables + offset);
2323 const MFloat DQ = vars[IP1] - vars[I];
2324 const MFloat DQP1 = vars[IP2] - vars[IP1];
2325 const MFloat DQM1 = vars[I] - vars[IM1];
2327 const MFloat ri = DQM1 / DQ;
2328 const MFloat rip = DQ / DQP1;
2333 m_QLeft[var] = vars[I] + (DQ * DSM1 + DQM1 * DS) * DSM * phii;
2334 m_QRight[var] = vars[IP1] - (DQP1 * DS + DQ * DSP1) * DSP * phiip;
2348 const MInt IM1 = I - IJK[dim];
2350 for(
MInt v = 0; v <
CV->noVariables; ++v) {
2376 for(
MInt dim = 0; dim <
nDim; dim++) {
2382 const MInt IP1 = I + IJK[dim];
2383 const MInt IM1 = I - IJK[dim];
2384 const MInt IP2 = I + 2 * IJK[dim];
2390 const MFloat DSP1 = sqrt(
POW2(x[IP2] - x[IP1]) +
POW2(
y[IP2] -
y[IP1]) +
POW2(z[IP2] - z[IP1]));
2394 const MFloat pIM2 = pvars[
PV->P][IM1];
2395 const MFloat pIM1 = pvars[
PV->P][I];
2396 const MFloat pIP2 = pvars[
PV->P][IP2];
2397 const MFloat pIP1 = pvars[
PV->P][IP1];
2399 const MFloat smps = DS * DSP1;
2400 const MFloat dummy = fabs(pIM2 - F2 * pIM1 + pIP1) / (pIM2 + F2 * pIM1 + pIP1);
2401 const MFloat dummy1 = fabs(pIM1 - F2 * pIP1 + pIP2) / (pIM1 + F2 * pIP1 + pIP2);
2405 for(
MInt var = 0; var <
PV->noVariables; ++var) {
2406 const MFloat DQ = pvars[var][IP1] - pvars[var][I];
2407 const MFloat DQP1 = pvars[var][IP2] - pvars[var][IP1];
2408 const MFloat DQM1 = pvars[var][I] - pvars[var][IM1];
2412 -
mMax(F0, (DQP1 * DQM1 * smps + F1B2 * epsLim) / (
POW2(DQP1 * DS) +
POW2(DQM1 * DSP1) + epsLim)))
2415 m_QLeft[var] = pvars[var][I] + DSM * (DSM1 * DQ + DS * DQM1) * phi;
2416 m_QRight[var] = pvars[var][IP1] - DSP * (DS * DQP1 + DSP1 * DQ) * phi;
2425 for(
MInt v = 0; v <
CV->noVariables; v++) {
2430 const MInt IM1 = I - IJK[dim];
2473 for(
MInt var = 0; var <
nDim + 2; var++) {
2474 const MUint offset = var * noCells;
2475 const MFloat*
const RESTRICT pvars = ALIGNED_F(cellVariables + offset);
2477 MFloat minNghbrDelta = F0;
2478 MFloat maxNghbrDelta = F0;
2479 MFloat effNghbrDelta = F0;
2482 for(
MInt dim1 = 0; dim1 <
nDim; dim1++) {
2483 for(
MInt rcnstructnNghbr = 0; rcnstructnNghbr < 2; rcnstructnNghbr++) {
2485 const MInt IP1 = cellId + IJK[dim1];
2486 const MInt IM1 = cellId - IJK[dim1];
2488 MFloat rcnstrctnNghbrValue = F0;
2489 if(rcnstructnNghbr == 0) {
2490 rcnstrctnNghbrValue = pvars[IP1];
2492 rcnstrctnNghbrValue = pvars[IM1];
2495 const MFloat tmpDelta = rcnstrctnNghbrValue - pvars[cellId];
2496 maxNghbrDelta =
mMax(maxNghbrDelta, tmpDelta);
2497 minNghbrDelta =
mMin(minNghbrDelta, tmpDelta);
2501 effNghbrDelta =
mMin(maxNghbrDelta, abs(minNghbrDelta));
2505 for(
MInt dim1 = 0; dim1 <
nDim; dim1++) {
2507 const MInt IP1 = cellId + IJK[dim1];
2508 const MInt IM1 = cellId - IJK[dim1];
2510 const MFloat DS = sqrt(
POW2(x[IP1] - x[cellId]) +
POW2(
y[IP1] -
y[cellId]) +
POW2(z[IP1] - z[cellId]));
2512 const MFloat DSM1 = sqrt(
POW2(x[cellId] - x[IM1]) +
POW2(
y[cellId] -
y[IM1]) +
POW2(z[cellId] - z[IM1]));
2518 * (DSM1 * sqrt(
POW2(x[IP1] - x[cellId]) +
POW2(
y[IP1] -
y[cellId]) +
POW2(z[IP1] - z[cellId]))
2519 + DS * sqrt(
POW2(x[cellId] - x[IM1]) +
POW2(
y[cellId] -
y[IM1]) +
POW2(z[cellId] - z[IM1])));
2521 (pvars[IP1] - pvars[IM1]) / sqrt(
POW2(x[IP1] - x[IM1]) +
POW2(
y[IP1] -
y[IM1]) +
POW2(z[IP1] - z[IM1]));
2522 srfcDelta += abs(DQ * dx1);
2536 for(
MInt dim1 = 0; dim1 <
nDim; dim1++) {
2537 for(
MInt rcnstructnNghbr = 0; rcnstructnNghbr < 2; rcnstructnNghbr++) {
2539 const MInt IP1 = cellId + IJK[dim1];
2540 const MInt IP2 = cellId + 2 * IJK[dim1];
2542 MFloat rcnstrctnNghbrValue = F0;
2543 if(rcnstructnNghbr == 0) {
2544 rcnstrctnNghbrValue = pvars[IP2];
2546 rcnstrctnNghbrValue = pvars[cellId];
2549 const MFloat tmpDelta = rcnstrctnNghbrValue - pvars[IP1];
2551 maxNghbrDelta =
mMax(maxNghbrDelta, tmpDelta);
2552 minNghbrDelta =
mMin(minNghbrDelta, tmpDelta);
2556 effNghbrDelta =
mMin(maxNghbrDelta, abs(minNghbrDelta));
2560 for(
MInt dim1 = 0; dim1 <
nDim; dim1++) {
2562 const MInt IP1 = cellId + IJK[dim1];
2563 const MInt IP2 = cellId + 2 * IJK[dim1];
2565 const MFloat DS = sqrt(
POW2(x[IP1] - x[cellId]) +
POW2(
y[IP1] -
y[cellId]) +
POW2(z[IP1] - z[cellId]));
2567 const MFloat DSP1 = sqrt(
POW2(x[IP2] - x[IP1]) +
POW2(
y[IP2] -
y[IP1]) +
POW2(z[IP2] - z[IP1]));
2574 * (DS * sqrt(
POW2(x[IP2] - x[IP1]) +
POW2(
y[IP2] -
y[IP1]) +
POW2(z[IP2] - z[IP1]))
2575 + DSP1 * sqrt(
POW2(x[IP1] - x[cellId]) +
POW2(
y[IP1] -
y[cellId]) +
POW2(z[IP1] - z[cellId])));
2576 const MFloat DQ = (pvars[IP2] - pvars[cellId])
2577 / sqrt(
POW2(x[IP2] - x[cellId]) +
POW2(
y[IP2] -
y[cellId]) +
POW2(z[IP2] - z[cellId]));
2579 srfcDelta += abs(DQ * dx2);
2586 for(
MInt dim1 = 0; dim1 <
nDim; dim1++) {
2588 const MInt IP1 = cellId + IJK[dim1];
2589 const MInt IM1 = cellId - IJK[dim1];
2590 const MInt IP2 = cellId + 2 * IJK[dim1];
2592 const MFloat DS = sqrt(
POW2(x[IP1] - x[cellId]) +
POW2(
y[IP1] -
y[cellId]) +
POW2(z[IP1] - z[cellId]));
2594 const MFloat DSM1 = sqrt(
POW2(x[cellId] - x[IM1]) +
POW2(
y[cellId] -
y[IM1]) +
POW2(z[cellId] - z[IM1]));
2595 const MFloat DSP1 = sqrt(
POW2(x[IP2] - x[IP1]) +
POW2(
y[IP2] -
y[IP1]) +
POW2(z[IP2] - z[IP1]));
2599 QLeft(cellId, var, dim1) =
2601 + DSM * (DSM1 * (pvars[IP1] - pvars[cellId]) + DS * (pvars[cellId] - pvars[IM1])) * minPhi(var, 0);
2602 QRight(cellId, var, dim1) =
2604 - DSP * (DS * (pvars[IP2] - pvars[IP1]) + DSP1 * (pvars[IP1] - pvars[cellId])) * minPhi(var, 1);
2612 for(
MInt dim = 0; dim <
nDim; dim++) {
2619 for(
MInt v = 0; v <
PV->noVariables; v++) {
2620 m_QLeft[v] = QLeft(cellId, v, dim);
2621 m_QRight[v] = QRight(cellId, v, dim);
2631 for(
MInt v = 0; v <
CV->noVariables; v++) {
2636 const MInt IM1 = I - IJK[dim];
2655 minPhi(var, cellPos) = (pow(effNghbrDelta, F2) + epsSqr + F2 * effNghbrDelta * srfcDelta)
2656 / (pow(effNghbrDelta, F2) + F2 * pow(srfcDelta, F2) + effNghbrDelta * srfcDelta + epsSqr);
2667 const MFloat eps = 1e-12;
2668 MFloat yps1 = effNghbrDelta / (srfcDelta + eps);
2669 minPhi(var, cellPos) =
mMin((yps1 * yps1 + F2 * yps1) / (yps1 * yps1 + yps1 + F2), F1);
2680 const MFloat eps = 1e-12;
2681 MFloat phi_max = effNghbrDelta / (srfcDelta + eps);
2682 minPhi(var, cellPos) =
mMin(phi_max, F1);
2696 const MFloat gammaMinusOne = gamma - 1.0;
2697 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
2710 const MFloat RHOL = QLeft[
PV->RHO];
2716 const MFloat RHOR = QRight[
PV->RHO];
2720 const MFloat fMetricLength = F1 / metricLength;
2723 const MFloat UUL = ((UL * surf0 + VL * surf1 + WL * surf2) - dxdtau) * fMetricLength;
2726 const MFloat UUR = ((UR * surf0 + VR * surf1 + WR * surf2) - dxdtau) * fMetricLength;
2733 const MFloat MAL = UUL / AL;
2734 const MFloat MAR = UUR / AR;
2736 const MFloat MALR = F1B2 * (MAL + MAR);
2739 const MFloat RHO_AL = RHOL * AL;
2740 const MFloat RHO_AR = RHOR * AR;
2742 const MFloat PLfRHOL = PL / RHOL;
2743 const MFloat PRfRHOR = PR / RHOR;
2745 const MFloat e0 = PLfRHOL * FgammaMinusOne + 0.5 * (
POW2(UL) +
POW2(VL) +
POW2(WL)) + PLfRHOL;
2746 const MFloat e1 = PRfRHOR * FgammaMinusOne + 0.5 * (
POW2(UR) +
POW2(VR) +
POW2(WR)) + PRfRHOR;
2748 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * metricLength;
2749 const MFloat RHOU2 = F1B2 * RHOU;
2751 const MFloat AbsRHO_U2 = fabs(RHOU2);
2755 flux[
CV->RHO_U][I] = RHOU2 * (UL + UR) + AbsRHO_U2 * (UL - UR) + PLR * surf0;
2756 flux[
CV->RHO_V][I] = RHOU2 * (VL + VR) + AbsRHO_U2 * (VL - VR) + PLR * surf1;
2757 flux[
CV->RHO_W][I] = RHOU2 * (WL + WR) + AbsRHO_U2 * (WL - WR) + PLR * surf2;
2758 flux[
CV->RHO_E][I] = RHOU2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1) + PLR * dxdtau;
2759 flux[
CV->RHO][I] = RHOU;
2788 const MFloat RHOL = QLeft[
PV->RHO];
2794 const MFloat RHOR = QRight[
PV->RHO];
2798 const MFloat fMetricLength = F1 / metricLength;
2801 const MFloat UUL = ((UL * surf0 + VL * surf1 + WL * surf2) - dxdtau) * fMetricLength;
2804 const MFloat UUR = ((UR * surf0 + VR * surf1 + WR * surf2) - dxdtau) * fMetricLength;
2810 const MFloat MAL = UUL / AL;
2811 const MFloat MAR = UUR / AR;
2813 const MFloat MALR = F1B2 * (MAL + MAR);
2831 const MFloat p4I4 = F4 * (p[IPJK] + p[IMJK]) - F6 * (p[I]) - p[IP2JK] - p[IM2JK];
2832 const MFloat p4J4 = F4 * (p[IJPK] + p[IJMK]) - F6 * (p[I]) - p[IJP2K] - p[IJM2K];
2833 const MFloat p4K4 = F4 * (p[IJKP] + p[IJKM]) - F6 * (p[I]) - p[IJKP2] - p[IJKM2];
2835 const MFloat pfac = fabs(p4I4) + fabs(p4J4) + fabs(p4K4);
2836 const MFloat facl = 1.0 / 1.3 * pfac;
2837 const MFloat fac = min(1 / 128.0, facl);
2839 const MFloat PLR = PL * (F1B2 + fac * MAL) + PR * (F1B2 - fac * MAR);
2841 const MFloat RHO_AL = RHOL * AL;
2842 const MFloat RHO_AR = RHOR * AR;
2844 const MFloat PLfRHOL = PL / RHOL;
2845 const MFloat PRfRHOR = PR / RHOR;
2847 const MFloat e0 = PLfRHOL * FgammaMinusOne + 0.5 * (
POW2(UL) +
POW2(VL) +
POW2(WL)) + PLfRHOL;
2848 const MFloat e1 = PRfRHOR * FgammaMinusOne + 0.5 * (
POW2(UR) +
POW2(VR) +
POW2(WR)) + PRfRHOR;
2850 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * metricLength;
2851 const MFloat RHOU2 = F1B2 * RHOU;
2853 const MFloat AbsRHO_U2 = fabs(RHOU2);
2856 flux[
CV->RHO_U][I] = RHOU2 * (UL + UR) + AbsRHO_U2 * (UL - UR) + PLR * surf0;
2857 flux[
CV->RHO_V][I] = RHOU2 * (VL + VR) + AbsRHO_U2 * (VL - VR) + PLR * surf1;
2858 flux[
CV->RHO_W][I] = RHOU2 * (WL + WR) + AbsRHO_U2 * (WL - WR) + PLR * surf2;
2859 flux[
CV->RHO_E][I] = RHOU2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1) + PLR * dxdtau;
2860 flux[
CV->RHO][I] = RHOU;
2865 const MFloat gammaMinusOne = gamma - 1.0;
2866 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
2875 const MFloat RHOL = QLeft[
PV->RHO];
2876 const MFloat FRHOL = F1 / RHOL;
2883 const MFloat RHOR = QRight[
PV->RHO];
2884 const MFloat FRHOR = F1 / RHOR;
2890 const MFloat PLfRHOL = PL / RHOL;
2891 const MFloat PRfRHOR = PR / RHOR;
2892 const MFloat e0 = PLfRHOL * FgammaMinusOne + 0.5 * (
POW2(UL) +
POW2(VL) +
POW2(WL)) + PLfRHOL;
2893 const MFloat e1 = PRfRHOR * FgammaMinusOne + 0.5 * (
POW2(UR) +
POW2(VR) +
POW2(WR)) + PRfRHOR;
2898 const MFloat FDGRAD = F1 / DGRAD;
2901 const MFloat UUL = ((UL * surf0 + VL * surf1 + WL * surf2) - dxdtau) * FDGRAD;
2904 const MFloat UUR = ((UR * surf0 + VR * surf1 + WR * surf2) - dxdtau) * FDGRAD;
2909 const MFloat FALR = 2.0 / (AL + AR);
2910 const MFloat ALPHAL = AL * FALR;
2911 const MFloat ALPHAR = AR * FALR;
2913 AL = sqrt(gamma * AL);
2914 AR = sqrt(gamma * AR);
2918 const MFloat XMAL = UUL / AL;
2919 const MFloat XMAR = UUR / AR;
2924 const MFloat RHOAL = AL * RHOL;
2925 const MFloat RHOAR = AR * RHOR;
2928 const MInt IP1 = I + IJK[dim];
2935 const MFloat SV1 = F0 * SV * DXDXEZ;
2936 const MFloat SV2 = F0 * SV * DYDXEZ;
2937 const MFloat SV3 = F0 * SV * DZDXEZ;
2942 MFloat FXMA = F1B2 * (XMAL1 + fabs(XMAL1));
2943 const MFloat XMALP = ALPHAL * (F1B4 *
POW2(XMAL1 + F1) - FXMA) + FXMA + (
mMax(F1, XMAL) - F1);
2944 FXMA = F1B2 * (XMAR1 - fabs(XMAR1));
2945 const MFloat XMARM = ALPHAR * (-F1B4 *
POW2(XMAR1 - F1) - FXMA) + FXMA + (
mMin(-F1, XMAR) + F1);
2947 const MFloat FLP = PL * ((F2 - XMAL1) *
POW2(F1 + XMAL1));
2948 const MFloat FRP = PR * ((F2 + XMAR1) *
POW2(F1 - XMAR1));
2949 const MFloat PLR = F1B4 * (FLP + FRP);
2951 const MFloat RHOUL = XMALP * RHOAL;
2952 const MFloat RHOUR = XMARM * RHOAR;
2953 const MFloat RHOU = RHOUL + RHOUR;
2954 const MFloat RHOU2 = F1B2 * RHOU;
2955 const MFloat ARHOU2 = fabs(RHOU2);
2957 const MFloat UUL2 = SV1 * UUL;
2958 const MFloat UUR2 = SV1 * UUR;
2961 const MFloat UUL3 = SV2 * UUL;
2962 const MFloat UUR3 = SV2 * UUR;
2965 const MFloat UUL4 = SV3 * UUL;
2966 const MFloat UUR4 = SV3 * UUR;
2972 flux[
CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0 + RHOUL * UUL2 + RHOUR * UUR2;
2973 flux[
CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1 + RHOUL * UUL3 + RHOUR * UUR3;
2974 flux[
CV->RHO_W][I] = RHOU2 * (WL + WR) + ARHOU2 * (WL - WR) + PLR * surf2 + RHOUL * UUL4 + RHOUR * UUR4;
2975 flux[
CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1) + PLR * dxdtau;
2976 flux[
CV->RHO][I] = RHOU;
2979template <MInt noVars>
2995 const MFloat gammaMinusOne = gamma - 1.0;
2996 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
3001 for(
MInt dim = 0; dim <
nDim; dim++) {
3004 const MInt IP1 = I + IJK[dim];
3005 dss[dim][I] = sqrt(
POW2(x[IP1] - x[I]) +
POW2(
y[IP1] -
y[I]) +
POW2(z[IP1] - z[I]));
3013 for(
MInt dim = 0; dim <
nDim; dim++) {
3016 const MInt IP1 = I + IJK[dim];
3017 const MInt IM1 = I - IJK[dim];
3018 const MInt IP2 = I + 2 * IJK[dim];
3020 const MFloat DS = dss[dim][I];
3021 const MFloat DSM1 = dss[dim][IM1];
3022 const MFloat DSP1 = dss[dim][IP1];
3029 const MFloat DQU = vars[
PV->U][IP1] - vars[
PV->U][I];
3030 const MFloat DQPU = vars[
PV->U][IP2] - vars[
PV->U][IP1];
3031 const MFloat DQMU = vars[
PV->U][I] - vars[
PV->U][IM1];
3032 const MFloat UL = vars[
PV->U][I] + DSM * (DSM1 * DQU + DS * DQMU);
3033 const MFloat UR = vars[
PV->U][IP1] - DSP * (DS * DQPU + DSP1 * DQU);
3035 const MFloat DQV = vars[
PV->V][IP1] - vars[
PV->V][I];
3036 const MFloat DQPV = vars[
PV->V][IP2] - vars[
PV->V][IP1];
3037 const MFloat DQMV = vars[
PV->V][I] - vars[
PV->V][IM1];
3038 const MFloat VL = vars[
PV->V][I] + DSM * (DSM1 * DQV + DS * DQMV);
3039 const MFloat VR = vars[
PV->V][IP1] - DSP * (DS * DQPV + DSP1 * DQV);
3041 const MFloat DQW = vars[
PV->W][IP1] - vars[
PV->W][I];
3042 const MFloat DQPW = vars[
PV->W][IP2] - vars[
PV->W][IP1];
3043 const MFloat DQMW = vars[
PV->W][I] - vars[
PV->W][IM1];
3044 const MFloat WL = vars[
PV->W][I] + DSM * (DSM1 * DQW + DS * DQMW);
3045 const MFloat WR = vars[
PV->W][IP1] - DSP * (DS * DQPW + DSP1 * DQW);
3047 const MFloat DQP = vars[
PV->P][IP1] - vars[
PV->P][I];
3048 const MFloat DQPP = vars[
PV->P][IP2] - vars[
PV->P][IP1];
3049 const MFloat DQMP = vars[
PV->P][I] - vars[
PV->P][IM1];
3050 const MFloat PL = vars[
PV->P][I] + DSM * (DSM1 * DQP + DS * DQMP);
3051 const MFloat PR = vars[
PV->P][IP1] - DSP * (DS * DQPP + DSP1 * DQP);
3053 const MFloat DQRHO = vars[
PV->RHO][IP1] - vars[
PV->RHO][I];
3054 const MFloat DQPRHO = vars[
PV->RHO][IP2] - vars[
PV->RHO][IP1];
3055 const MFloat DQMRHO = vars[
PV->RHO][I] - vars[
PV->RHO][IM1];
3056 const MFloat RHOL = vars[
PV->RHO][I] + DSM * (DSM1 * DQRHO + DS * DQMRHO);
3057 const MFloat RHOR = vars[
PV->RHO][IP1] - DSP * (DS * DQPRHO + DSP1 * DQRHO);
3066 const MFloat fMetricLength = F1 / metricLength;
3069 const MFloat UUL = ((UL * surf0 + VL * surf1 + WL * surf2) - dxdtau) * fMetricLength;
3072 const MFloat UUR = ((UR * surf0 + VR * surf1 + WR * surf2) - dxdtau) * fMetricLength;
3079 const MFloat MAL = UUL / AL;
3080 const MFloat MAR = UUR / AR;
3082 const MFloat MALR = F1B2 * (MAL + MAR);
3085 const MFloat RHO_AL = RHOL * AL;
3086 const MFloat RHO_AR = RHOR * AR;
3088 const MFloat PLfRHOL = PL / RHOL;
3089 const MFloat PRfRHOR = PR / RHOR;
3091 const MFloat e0 = PLfRHOL * FgammaMinusOne + 0.5 * (
POW2(UL) +
POW2(VL) +
POW2(WL)) + PLfRHOL;
3092 const MFloat e1 = PRfRHOR * FgammaMinusOne + 0.5 * (
POW2(UR) +
POW2(VR) +
POW2(WR)) + PRfRHOR;
3094 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * metricLength;
3095 const MFloat RHOU2 = F1B2 * RHOU;
3097 const MFloat AbsRHO_U2 = fabs(RHOU2);
3099 flux[
CV->RHO_U][I] = RHOU2 * (UL + UR) + AbsRHO_U2 * (UL - UR) + PLR * surf0;
3100 flux[
CV->RHO_V][I] = RHOU2 * (VL + VR) + AbsRHO_U2 * (VL - VR) + PLR * surf1;
3101 flux[
CV->RHO_W][I] = RHOU2 * (WL + WR) + AbsRHO_U2 * (WL - WR) + PLR * surf2;
3102 flux[
CV->RHO_E][I] = RHOU2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1) + PLR * dxdtau;
3103 flux[
CV->RHO][I] = RHOU;
3107 for(
MUint v = 0; v < noVars; v++) {
3110 const MInt IM1 = I - IJK[dim];
3111 const MUint offset = v * noCells;
3112 MFloat*
const RESTRICT
rhs = ALIGNED_F(cellRhs + offset);
3113 rhs[I] += flux[v][IM1] - flux[v][I];
3119template void FvStructuredSolver3D::Muscl_AusmLES<5>();
3120template void FvStructuredSolver3D::Muscl_AusmLES<6>();
3121template void FvStructuredSolver3D::Muscl_AusmLES<7>();
3123template <MInt noVars>
3141 const MFloat gammaMinusOne = gamma - 1.0;
3142 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
3153 for(
MInt dim = 0; dim <
nDim; dim++) {
3155#pragma omp parallel for
3157 for(
MInt k = 0; k < noCellsKP1; k++) {
3158 for(
MInt j = 0; j < noCellsJP1; j++) {
3159 for(
MInt i = 0; i < noCellsIP1; i++) {
3161 const MInt IP1 = I + IJK[dim];
3162 dss[dim][I] = sqrt(
POW2(x[IP1] - x[I]) +
POW2(
y[IP1] -
y[I]) +
POW2(z[IP1] - z[I]));
3171 for(
MInt dim = 0; dim <
nDim; dim++) {
3173#pragma omp parallel for
3175 for(
MInt k = 1; k < noCellsK; k++) {
3176 for(
MInt j = 1; j < noCellsJ; j++) {
3177#if defined(MAIA_INTEL_COMPILER)
3179#pragma vector always
3181 for(
MInt i = 1; i < noCellsI; i++) {
3183 const MInt IP1 = I + IJK[dim];
3184 const MInt IM1 = I - IJK[dim];
3185 const MInt IP2 = I + 2 * IJK[dim];
3187 const MFloat DS = dss[dim][I];
3188 const MFloat DSM1 = dss[dim][IM1];
3189 const MFloat DSP1 = dss[dim][IP1];
3196 const MFloat DQU = vars[
PV->U][IP1] - vars[
PV->U][I];
3197 const MFloat DQPU = vars[
PV->U][IP2] - vars[
PV->U][IP1];
3198 const MFloat DQMU = vars[
PV->U][I] - vars[
PV->U][IM1];
3199 const MFloat UL = vars[
PV->U][I] + DSM * (DSM1 * DQU + DS * DQMU);
3200 const MFloat UR = vars[
PV->U][IP1] - DSP * (DS * DQPU + DSP1 * DQU);
3202 const MFloat DQV = vars[
PV->V][IP1] - vars[
PV->V][I];
3203 const MFloat DQPV = vars[
PV->V][IP2] - vars[
PV->V][IP1];
3204 const MFloat DQMV = vars[
PV->V][I] - vars[
PV->V][IM1];
3205 const MFloat VL = vars[
PV->V][I] + DSM * (DSM1 * DQV + DS * DQMV);
3206 const MFloat VR = vars[
PV->V][IP1] - DSP * (DS * DQPV + DSP1 * DQV);
3208 const MFloat DQW = vars[
PV->W][IP1] - vars[
PV->W][I];
3209 const MFloat DQPW = vars[
PV->W][IP2] - vars[
PV->W][IP1];
3210 const MFloat DQMW = vars[
PV->W][I] - vars[
PV->W][IM1];
3211 const MFloat WL = vars[
PV->W][I] + DSM * (DSM1 * DQW + DS * DQMW);
3212 const MFloat WR = vars[
PV->W][IP1] - DSP * (DS * DQPW + DSP1 * DQW);
3214 const MFloat DQP = vars[
PV->P][IP1] - vars[
PV->P][I];
3215 const MFloat DQPP = vars[
PV->P][IP2] - vars[
PV->P][IP1];
3216 const MFloat DQMP = vars[
PV->P][I] - vars[
PV->P][IM1];
3217 const MFloat PL = vars[
PV->P][I] + DSM * (DSM1 * DQP + DS * DQMP);
3218 const MFloat PR = vars[
PV->P][IP1] - DSP * (DS * DQPP + DSP1 * DQP);
3220 const MFloat DQRHO = vars[
PV->RHO][IP1] - vars[
PV->RHO][I];
3221 const MFloat DQPRHO = vars[
PV->RHO][IP2] - vars[
PV->RHO][IP1];
3222 const MFloat DQMRHO = vars[
PV->RHO][I] - vars[
PV->RHO][IM1];
3223 const MFloat RHOL = vars[
PV->RHO][I] + DSM * (DSM1 * DQRHO + DS * DQMRHO);
3224 const MFloat RHOR = vars[
PV->RHO][IP1] - DSP * (DS * DQPRHO + DSP1 * DQRHO);
3233 const MFloat fMetricLength = F1 / metricLength;
3236 const MFloat UUL = ((UL * surf0 + VL * surf1 + WL * surf2) - dxdtau) * fMetricLength;
3239 const MFloat UUR = ((UR * surf0 + VR * surf1 + WR * surf2) - dxdtau) * fMetricLength;
3246 const MFloat MAL = UUL / AL;
3247 const MFloat MAR = UUR / AR;
3249 const MFloat MALR = F1B2 * (MAL + MAR);
3267 const MFloat p4I4 = F4 * (p[IPJK] + p[IMJK]) - F6 * (p[I]) - p[IP2JK] - p[IM2JK];
3268 const MFloat p4J4 = F4 * (p[IJPK] + p[IJMK]) - F6 * (p[I]) - p[IJP2K] - p[IJM2K];
3269 const MFloat p4K4 = F4 * (p[IJKP] + p[IJKM]) - F6 * (p[I]) - p[IJKP2] - p[IJKM2];
3271 const MFloat pfac = fabs(p4I4) + fabs(p4J4) + fabs(p4K4);
3272 const MFloat facl = 1.0 / 1.3 * pfac;
3273 const MFloat fac = min(1.0 / 128.0, facl);
3275 const MFloat PLR = PL * (F1B2 + fac * MAL) + PR * (F1B2 - fac * MAR);
3277 const MFloat RHO_AL = RHOL * AL;
3278 const MFloat RHO_AR = RHOR * AR;
3280 const MFloat PLfRHOL = PL / RHOL;
3281 const MFloat PRfRHOR = PR / RHOR;
3283 const MFloat e0 = PLfRHOL * FgammaMinusOne + 0.5 * (
POW2(UL) +
POW2(VL) +
POW2(WL)) + PLfRHOL;
3284 const MFloat e1 = PRfRHOR * FgammaMinusOne + 0.5 * (
POW2(UR) +
POW2(VR) +
POW2(WR)) + PRfRHOR;
3286 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * metricLength;
3287 const MFloat RHOU2 = F1B2 * RHOU;
3289 const MFloat AbsRHO_U2 = fabs(RHOU2);
3291 flux[
CV->RHO_U][I] = RHOU2 * (UL + UR) + AbsRHO_U2 * (UL - UR) + PLR * surf0;
3292 flux[
CV->RHO_V][I] = RHOU2 * (VL + VR) + AbsRHO_U2 * (VL - VR) + PLR * surf1;
3293 flux[
CV->RHO_W][I] = RHOU2 * (WL + WR) + AbsRHO_U2 * (WL - WR) + PLR * surf2;
3294 flux[
CV->RHO_E][I] = RHOU2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1) + PLR * dxdtau;
3295 flux[
CV->RHO][I] = RHOU;
3301 for(
MUint v = 0; v < noVars; v++) {
3303#pragma omp parallel for
3307#if defined(MAIA_INTEL_COMPILER)
3309#pragma vector always
3313 const MInt IM1 = I - IJK[dim];
3314 const MUint offset = v * noCells;
3315 MFloat*
const RESTRICT
rhs = ALIGNED_F(cellRhs + offset);
3316 rhs[I] += flux[v][IM1] - flux[v][I];
3324template void FvStructuredSolver3D::Muscl_AusmLES_PTHRC<5>();
3325template void FvStructuredSolver3D::Muscl_AusmLES_PTHRC<6>();
3326template void FvStructuredSolver3D::Muscl_AusmLES_PTHRC<7>();
3329template <MInt noVars>
3346 const MFloat gammaMinusOne = gamma - 1.0;
3347 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
3358 for(
MInt dim = 0; dim <
nDim; dim++) {
3360#pragma omp parallel for
3362 for(
MInt k = 0; k < noCellsKP1; k++) {
3363 for(
MInt j = 0; j < noCellsJP1; j++) {
3364 for(
MInt i = 0; i < noCellsIP1; i++) {
3366 const MInt IP1 = I + IJK[dim];
3367 dss[dim][I] = sqrt(
POW2(x[IP1] - x[I]) +
POW2(
y[IP1] -
y[I]) +
POW2(z[IP1] - z[I]));
3377 for(
MInt dim = 0; dim <
nDim; dim++) {
3379#pragma omp parallel for
3381 for(
MInt k = 1; k < noCellsK; k++) {
3382 for(
MInt j = 1; j < noCellsJ; j++) {
3383#if defined(MAIA_INTEL_COMPILER)
3385#pragma vector always
3387 for(
MInt i = 1; i < noCellsI; i++) {
3389 const MInt IP1 = I + IJK[dim];
3390 const MInt IM1 = I - IJK[dim];
3391 const MInt IP2 = I + 2 * IJK[dim];
3393 const MFloat DS = dss[dim][I];
3394 const MFloat DSM1 = dss[dim][IM1];
3395 const MFloat DSP1 = dss[dim][IP1];
3402 const MFloat DQU = vars[
PV->U][IP1] - vars[
PV->U][I];
3403 const MFloat DQPU = vars[
PV->U][IP2] - vars[
PV->U][IP1];
3404 const MFloat DQMU = vars[
PV->U][I] - vars[
PV->U][IM1];
3405 MFloat UL = vars[
PV->U][I] + DSM * (DSM1 * DQU + DS * DQMU);
3406 MFloat UR = vars[
PV->U][IP1] - DSP * (DS * DQPU + DSP1 * DQU);
3408 const MFloat DQV = vars[
PV->V][IP1] - vars[
PV->V][I];
3409 const MFloat DQPV = vars[
PV->V][IP2] - vars[
PV->V][IP1];
3410 const MFloat DQMV = vars[
PV->V][I] - vars[
PV->V][IM1];
3411 MFloat VL = vars[
PV->V][I] + DSM * (DSM1 * DQV + DS * DQMV);
3412 MFloat VR = vars[
PV->V][IP1] - DSP * (DS * DQPV + DSP1 * DQV);
3414 const MFloat DQW = vars[
PV->W][IP1] - vars[
PV->W][I];
3415 const MFloat DQPW = vars[
PV->W][IP2] - vars[
PV->W][IP1];
3416 const MFloat DQMW = vars[
PV->W][I] - vars[
PV->W][IM1];
3417 MFloat WL = vars[
PV->W][I] + DSM * (DSM1 * DQW + DS * DQMW);
3418 MFloat WR = vars[
PV->W][IP1] - DSP * (DS * DQPW + DSP1 * DQW);
3420 const MFloat DQP = vars[
PV->P][IP1] - vars[
PV->P][I];
3421 const MFloat DQPP = vars[
PV->P][IP2] - vars[
PV->P][IP1];
3422 const MFloat DQMP = vars[
PV->P][I] - vars[
PV->P][IM1];
3423 const MFloat PL = vars[
PV->P][I] + DSM * (DSM1 * DQP + DS * DQMP);
3424 const MFloat PR = vars[
PV->P][IP1] - DSP * (DS * DQPP + DSP1 * DQP);
3426 const MFloat DQRHO = vars[
PV->RHO][IP1] - vars[
PV->RHO][I];
3427 const MFloat DQPRHO = vars[
PV->RHO][IP2] - vars[
PV->RHO][IP1];
3428 const MFloat DQMRHO = vars[
PV->RHO][I] - vars[
PV->RHO][IM1];
3429 const MFloat RHOL = vars[
PV->RHO][I] + DSM * (DSM1 * DQRHO + DS * DQMRHO);
3430 const MFloat RHOR = vars[
PV->RHO][IP1] - DSP * (DS * DQPRHO + DSP1 * DQRHO);
3437 const MFloat FRHOL = F1 / RHOL;
3438 const MFloat FRHOR = F1 / RHOR;
3440 const MFloat PLfRHOL = PL / RHOL;
3441 const MFloat PRfRHOR = PR / RHOR;
3442 const MFloat e0 = PLfRHOL * FgammaMinusOne + 0.5 * (
POW2(UL) +
POW2(VL) +
POW2(WL)) + PLfRHOL;
3443 const MFloat e1 = PRfRHOR * FgammaMinusOne + 0.5 * (
POW2(UR) +
POW2(VR) +
POW2(WR)) + PRfRHOR;
3448 const MFloat FDGRAD = F1 / DGRAD;
3451 const MFloat UUL = ((UL * surf0 + VL * surf1 + WL * surf2) - dxdtau) * FDGRAD;
3454 const MFloat UUR = ((UR * surf0 + VR * surf1 + WR * surf2) - dxdtau) * FDGRAD;
3459 const MFloat FALR = 2.0 / (AL + AR);
3460 const MFloat ALPHAL = AL * FALR;
3461 const MFloat ALPHAR = AR * FALR;
3463 AL = sqrt(gamma * AL);
3464 AR = sqrt(gamma * AR);
3468 const MFloat XMAL = UUL / AL;
3469 const MFloat XMAR = UUR / AR;
3474 const MFloat RHOAL = AL * RHOL;
3475 const MFloat RHOAR = AR * RHOR;
3482 const MFloat SV1 = F0 * SV * DXDXEZ;
3483 const MFloat SV2 = F0 * SV * DYDXEZ;
3484 const MFloat SV3 = F0 * SV * DZDXEZ;
3489 MFloat FXMA = F1B2 * (XMAL1 + fabs(XMAL1));
3490 const MFloat XMALP = ALPHAL * (F1B4 *
POW2(XMAL1 + F1) - FXMA) + FXMA + (
mMax(F1, XMAL) - F1);
3491 FXMA = F1B2 * (XMAR1 - fabs(XMAR1));
3492 const MFloat XMARM = ALPHAR * (-F1B4 *
POW2(XMAR1 - F1) - FXMA) + FXMA + (
mMin(-F1, XMAR) + F1);
3494 const MFloat FLP = PL * ((F2 - XMAL1) *
POW2(F1 + XMAL1));
3495 const MFloat FRP = PR * ((F2 + XMAR1) *
POW2(F1 - XMAR1));
3496 const MFloat PLR = F1B4 * (FLP + FRP);
3498 const MFloat RHOUL = XMALP * RHOAL;
3499 const MFloat RHOUR = XMARM * RHOAR;
3500 const MFloat RHOU = RHOUL + RHOUR;
3501 const MFloat RHOU2 = F1B2 * RHOU;
3502 const MFloat ARHOU2 = fabs(RHOU2);
3504 const MFloat UUL2 = SV1 * UUL;
3505 const MFloat UUR2 = SV1 * UUR;
3508 const MFloat UUL3 = SV2 * UUL;
3509 const MFloat UUR3 = SV2 * UUR;
3512 const MFloat UUL4 = SV3 * UUL;
3513 const MFloat UUR4 = SV3 * UUR;
3517 flux[
CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0 + RHOUL * UUL2 + RHOUR * UUR2;
3518 flux[
CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1 + RHOUL * UUL3 + RHOUR * UUR3;
3519 flux[
CV->RHO_W][I] = RHOU2 * (WL + WR) + ARHOU2 * (WL - WR) + PLR * surf2 + RHOUL * UUL4 + RHOUR * UUR4;
3520 flux[
CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1) + PLR * dxdtau;
3521 flux[
CV->RHO][I] = RHOU;
3527 for(
MUint v = 0; v < noVars; v++) {
3529#pragma omp parallel for
3533#if defined(MAIA_INTEL_COMPILER)
3535#pragma vector always
3539 const MInt IM1 = I - IJK[dim];
3540 const MUint offset = v * noCells;
3541 MFloat*
const RESTRICT
rhs = ALIGNED_F(cellRhs + offset);
3542 rhs[I] += flux[v][IM1] - flux[v][I];
3550template void FvStructuredSolver3D::Muscl_AusmDV<5>();
3551template void FvStructuredSolver3D::Muscl_AusmDV<6>();
3552template void FvStructuredSolver3D::Muscl_AusmDV<7>();
3555template <FvStructuredSolver3D::fluxmethod ausm, MInt noVars>
3573 for(
MInt dim = 0; dim <
nDim; dim++) {
3575 const MFloat*
const RESTRICT length = ALIGNED_F(cellLength + dimOffset);
3581 const MInt IP1 = I + IJK[dim];
3582 const MInt IM1 = I - IJK[dim];
3583 const MInt IP2 = I + 2 * IJK[dim];
3585 const MFloat rp = (length[I] + length[IP1]) / (F2 * length[I]);
3586 const MFloat rm = (length[I] + length[IM1]) / (F2 * length[I]);
3587 const MFloat f = phi / (F2 * (rp + rm));
3588 const MFloat f1 = (rm + kappa * phi) / rp;
3589 const MFloat f2 = (rp - kappa * phi) / rm;
3591 const MFloat rp1 = (length[IP1] + length[IP2]) / (F2 * length[IP1]);
3592 const MFloat rm1 = (length[IP1] + length[I]) / (F2 * length[IP1]);
3593 const MFloat fa = phi / (F2 * (rp1 + rm1));
3594 const MFloat fb = (rm1 - kappa * phi) / rp1;
3595 const MFloat fc = (rp1 + kappa * phi) / rm1;
3597 for(
MUint v = 0; v < noVars; v++) {
3599 const MFloat*
const RESTRICT vars = ALIGNED_F(cellVariables + offset);
3601 const MFloat DQ = (vars[IP1] - vars[I]);
3602 const MFloat DQM1 = (vars[I] - vars[IM1]);
3603 qleft[v] = vars[I] + f * (f1 * DQ + f2 * DQM1);
3606 const MFloat DQP1 = (vars[IP2] - vars[IP1]);
3607 const MFloat DQ1 = (vars[IP1] - vars[I]);
3608 qright[v] = vars[IP1] - fa * (fb * DQP1 + fc * DQ1);
3617 for(
MUint v = 0; v < noVars; v++) {
3622 const MInt IM1 = I - IJK[dim];
3623 const MUint offset = v * noCells;
3624 MFloat*
const RESTRICT
rhs = ALIGNED_F(cellRhs + offset);
3625 rhs[I] += flux[v][IM1] - flux[v][I];
3633template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmLES, 5>();
3634template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmLES, 6>();
3635template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmLES, 7>();
3637template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmLES_PTHRC, 5>();
3638template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmLES_PTHRC, 6>();
3639template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmLES_PTHRC, 7>();
3641template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmDV, 5>();
3642template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmDV, 6>();
3643template void FvStructuredSolver3D::MusclStretched_<&FvStructuredSolver3D::AusmDV, 7>();
3655 m_grid->saveCellJacobian();
3699 MFloat minCoordinate = F0, maxCoordinate = F0, minCoordinateGlobal = F0, maxCoordinateGlobal = F0;
3700 MInt lowPoint = -1, highPoint = -1;
3716 "minCoordinateGlobal");
3718 "maxCoordinateGlobal");
3783 tripX = Context::getSolverProperty<MFloat>(
"tripXLength",
m_solverId, AT_, &tripX);
3817 tripY = Context::getSolverProperty<MFloat>(
"tripYHeight",
m_solverId, AT_, &tripY);
3836 tripZ = Context::getSolverProperty<MFloat>(
"tripCutoffZ",
m_solverId, AT_, &tripZ);
3888 timeCutoff = Context::getSolverProperty<MFloat>(
"tripDeltaTime",
m_solverId, AT_, &timeCutoff);
3895 mTerm(1,
"Properties RNGSeed or seedRNGWithTime not compatible with SandpaperTrip!");
3910 "localDomainWidth",
"m_tripDomainWidth");
3923 const MInt blockId = 0;
3932 localCoords.fill(-99999.0);
3933 localNormalVec.fill(-99999.0);
3942 const MInt localPointId = i;
3943 const MInt globalPointId = offset + i;
3949 MFloat normalVec[3] = {F0, F0, F0};
3950 for(
MInt dim = 0; dim <
nDim; dim++) {
3951 normalVec[dim] =
m_grid->m_coordinates[dim][pIJPK] -
m_grid->m_coordinates[dim][pIJK];
3955 const MFloat normalLength = sqrt(
POW2(normalVec[0]) +
POW2(normalVec[1]) +
POW2(normalVec[2]));
3956 for(
MInt dim = 0; dim <
nDim; dim++) {
3957 normalVec[dim] /= normalLength;
3960 localCoords(globalPointId) =
m_grid->m_coordinates[0][localPointId];
3964 localNormalVec(globalPointId) = normalVec[0];
3978 vector<Point<3>> pts;
3980 const MInt gIJK = i;
3981 const MInt gIPJK = i + 1;
3999 (void)tree.nearest(pt, distance);
4008 m_log <<
"=================================================" << endl
4009 <<
" SANDPAPER TRIP PROPERTIES " << endl
4010 <<
"=================================================" << endl;
4012 m_log <<
"######### TRIP NUMBER " << i <<
" ###########" << endl
4022 <<
"###########################################" << endl;
4024 m_log <<
"=================================================" << endl;
4040 stringstream tripPath;
4041 tripPath <<
"/trip";
4042 ParallelIo::size_type dummyOffset = 0;
4046 stringstream restartFileName;
4050 pio.readArray(
m_tripModesG, tripPath.str(),
"tripModesG", 1, &dummyOffset, &dataSize);
4051 pio.readArray(
m_tripModesH1, tripPath.str(),
"tripModesH1", 1, &dummyOffset, &dataSize);
4052 pio.readArray(
m_tripModesH2, tripPath.str(),
"tripModesH2", 1, &dummyOffset, &dataSize);
4084 const MFloat b = 3 * pow(p, 2) - 2 * pow(p, 3);
4107 const MFloat forceStrength =
4142 const MFloat b = 3 * pow(p, 2) - 2 * pow(p, 3);
4165 const MFloat forceStrength =
4180 if(ii == 0 &&
y > 0.0) {
4182 }
else if(ii == 1 &&
y < 0.0) {
4186 force = forceStrength * factor;
4190 const MInt globalPointId = offsetPoints + i;
4197 const MFloat u_force = wallNormalVec[0] * force;
4203 const MFloat v_force = wallNormalVec[1] * force;
4216 MFloat maxLocalValue = -999999.0;
4217 MFloat minLocalValue = 999999.0;
4219 MFloat* phik = &modes[noModes];
4221 for(
MInt k = 0; k < noCells; k++) {
4222 const MFloat z = coords[k];
4223 for(
MInt n = 0; n < noModes; n++) {
4224 forceCoef[k] += sin(z * ak[n] + phik[n]);
4227 maxLocalValue =
mMax(maxLocalValue, forceCoef[k]);
4228 minLocalValue =
mMin(minLocalValue, forceCoef[k]);
4232 MFloat maxGlobalValue = 0.0;
4233 MFloat minGlobalValue = 0.0;
4240 for(
MInt k = 0; k < noCells; k++) {
4241 forceCoef[k] = 2 * (forceCoef[k] - minGlobalValue) / (maxGlobalValue - minGlobalValue) - 1.0;
4248 const MFloat minWavenumber = 2 * PI / maxWaveLength;
4249 const MFloat maxWavenumber = 2 * PI / minWaveLength;
4252 MFloat* phik = &modes[noModes];
4254 for(
MInt n = 0; n < noModes; n++) {
4255 ak[n] = (maxWavenumber - minWavenumber) * (rand() /
MFloat(RAND_MAX)) + minWavenumber;
4256 phik[n] = 2 * PI * rand() /
MFloat(RAND_MAX);
4304 if(!variables.empty()) {
4305 MInt totalNoCellsIJ = (
m_grid->getMyBlockNoCells(2) *
m_grid->getMyBlockNoCells(1));
4306 const MInt noVars = variables.size();
4307 MFloatScratchSpace localSpannwiseVars(totalNoCellsIJ, noVars, AT_,
"localSpannwiseVars");
4308 MFloatScratchSpace globalSpannwiseVars(totalNoCellsIJ, noVars, AT_,
"globalSpannwiseVars");
4310 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4314 MInt localCellId3D =
4317 localSpannwiseVars(globalCellId2D, varPos) += variables[varPos][localCellId3D];
4323 MPI_Allreduce(&localSpannwiseVars(0, 0), &globalSpannwiseVars(0, 0), totalNoCellsIJ * noVars, MPI_DOUBLE, MPI_SUM,
4326 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4327 for(
MInt cellId = 0; cellId < totalNoCellsIJ; cellId++) {
4328 globalSpannwiseVars(cellId, varPos) /=
m_grid->getMyBlockNoCells(0);
4332 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4336 MInt localCellId3D =
4339 variables[varPos][localCellId3D] = globalSpannwiseVars(globalCellId2D, varPos);
4369 const MFloat lenXi = sqrt(
POW2(metric[0][cellId]) +
POW2(metric[1][cellId]) +
POW2(metric[2][cellId]));
4370 const MFloat lenEt = sqrt(
POW2(metric[3][cellId]) +
POW2(metric[4][cellId]) +
POW2(metric[5][cellId]));
4371 const MFloat lenZe = sqrt(
POW2(metric[6][cellId]) +
POW2(metric[7][cellId]) +
POW2(metric[8][cellId]));
4385 U_c -= dxtx[cellId];
4386 V_c -= dxty[cellId];
4387 W_c -= dxtz[cellId];
4394 const MFloat eigenvalue = U_c + V_c + W_c + speedOfSound * (lenXi + lenEt + lenZe);
4420 const MInt noVars =
CV->noVariables;
4433 for(
MInt v = 0; v < noVars; v++) {
4434 const MUint offset = v * noCells;
4435 MFloat*
const RESTRICT oldCellVars = ALIGNED_F(oldVars + offset);
4436 const MFloat*
const RESTRICT cellVars = ALIGNED_F(vars + offset);
4437 maia::parallelFor<true>(0,
m_noCells, [=](
MInt cellId) { oldCellVars[cellId] = cellVars[cellId]; });
4445 for(
MInt v = 0; v < noVars; v++) {
4446 const MUint cellOffset = v * noCells;
4447 MFloat*
const RESTRICT cellVars = vars + cellOffset;
4448 const MFloat*
const RESTRICT oldCellVars = oldVars + cellOffset;
4449 const MFloat*
const RESTRICT cellRhs =
rhs + cellOffset;
4455 cellVars[cellId] = oldCellVars[cellId] + factor * cellRhs[cellId];
4459 for(
MInt v = 0; v < noVars; v++) {
4460 const MUint cellOffset = v * noCells;
4461 MFloat*
const RESTRICT cellVars = vars + cellOffset;
4462 const MFloat*
const RESTRICT oldCellVars = oldVars + cellOffset;
4463 const MFloat*
const RESTRICT cellRhs =
rhs + cellOffset;
4468 (oldCellVars[cellId] * oldCellJac[cellId] + rkFactor * cellRhs[cellId]) / cellJac[cellId];
4472 for(
MInt v = 0; v < noVars; v++) {
4473 const MUint cellOffset = v * noCells;
4474 MFloat*
const RESTRICT cellVars = vars + cellOffset;
4475 const MFloat*
const RESTRICT oldCellVars = oldVars + cellOffset;
4476 const MFloat*
const RESTRICT cellRhs =
rhs + cellOffset;
4481 cellVars[cellId] = oldCellVars[cellId] + factor * cellRhs[cellId];
4488 for(
MInt v = 0; v < noVars; v++) {
4489 const MUint cellOffset = v * noCells;
4490 MFloat*
const RESTRICT cellVars = vars + cellOffset;
4491 const MFloat*
const RESTRICT oldCellVars = oldVars + cellOffset;
4492 const MFloat*
const RESTRICT cellRhs =
rhs + cellOffset;
4498 rkAlpha * cellVars[cellId] + (F1 - rkAlpha) * oldCellVars[cellId] - factor * cellRhs[cellId];
4504 stringstream errorMessage;
4505 errorMessage <<
"Given RungeKutta Order " <<
m_rungeKuttaOrder <<
" not implemented! " << endl;
4506 mTerm(1, AT_, errorMessage.str());
4526 cout <<
"enterin addDisturbance " << endl;
4527 cout <<
" m_time = " <<
m_time << endl;
4528 cout <<
"m_timeStep = " <<
m_timeStep << endl;
4531 MFloat pi2 = 8.0 * atan(1);
4540 cout <<
"sin " << amp1 << endl;
4541 MFloat centerloc = 3.1415;
4555 MFloat disturb = amp1 * exp(-(a1 + a2 + a3) /
POW2(r) / 2.0);
4566 cout <<
"leaving addDisturbance " << endl;
4599 m_grid->extrapolateGhostPointCoordinates();
4606 m_grid->computeCellCenterCoordinates();
4617 m_grid->writePartitionedGrid();
4636 MInt pos[3], fix[3], mirror[3], ijk[3], extendijk[3];
4637 MInt pointId, FixPointId, MirrorPointId;
4642 extendijk[index] = 0;
4644 for(ijk[2] = start[2]; ijk[2] < end[2] + extendijk[2]; ++ijk[2]) {
4645 for(ijk[1] = start[1]; ijk[1] < end[1] + extendijk[1]; ++ijk[1]) {
4646 for(ijk[0] = start[0]; ijk[0] < end[0] + extendijk[0]; ++ijk[0]) {
4647 for(
MInt m = 0; m < 3; ++m) {
4650 pos[m] = ijk[m] + 1;
4652 mirror[m] = 2 * fix[m] - pos[m];
4656 mirror[m] = 2 * fix[m] - pos[m];
4665 pointId =
pointIndex(pos[0], pos[1], pos[2]);
4666 FixPointId =
pointIndex(fix[0], fix[1], fix[2]);
4667 MirrorPointId =
pointIndex(mirror[0], mirror[1], mirror[2]);
4669 for(
MInt dim = 0; dim <
nDim; dim++) {
4670 m_grid->m_coordinates[dim][pointId] =
4671 (2 *
m_grid->m_coordinates[dim][FixPointId] -
m_grid->m_coordinates[dim][MirrorPointId]);
4686 for(
auto& snd : sndComm) {
4690 std::array<MInt, nDim> begin{snd->startInfoCells[0], snd->startInfoCells[1], snd->startInfoCells[2]};
4691 std::array<MInt, nDim> end{snd->endInfoCells[0], snd->endInfoCells[1], snd->endInfoCells[2]};
4692 std::array<MInt, nDim> size{snd->endInfoCells[0] - snd->startInfoCells[0],
4693 snd->endInfoCells[1] - snd->startInfoCells[1],
4694 snd->endInfoCells[2] - snd->startInfoCells[2]};
4695 const MInt totalSize = size[0] * size[1] * size[2];
4698 auto* cellBuffer = snd->cellBuffer.get();
4699 auto* variables = &(snd->variables[0]);
4700 for(
MInt var = 0; var < snd->noVars; var++) {
4701 maia::parallelFor<true, nDim>(begin, end, [=](
const MInt& i,
const MInt& j,
const MInt& k) {
4703 const MInt bufferId = totalSize * var + (i - begin[0]) + ((j - begin[1]) + (k - begin[2]) * size[1]) * size[0];
4704 cellBuffer[bufferId] = variables[var][cellId];
4717 for(
auto& rcv : rcvComm) {
4721 std::array<MInt, nDim> begin{rcv->startInfoCells[0], rcv->startInfoCells[1], rcv->startInfoCells[2]};
4722 std::array<MInt, nDim> end{rcv->endInfoCells[0], rcv->endInfoCells[1], rcv->endInfoCells[2]};
4723 std::array<MInt, nDim> size{rcv->endInfoCells[0] - rcv->startInfoCells[0],
4724 rcv->endInfoCells[1] - rcv->startInfoCells[1],
4725 rcv->endInfoCells[2] - rcv->startInfoCells[2]};
4726 const MInt totalSize = size[0] * size[1] * size[2];
4728 std::array<MInt, nDim> stepBuffer{0};
4729 std::array<MInt, nDim> startBuffer{0};
4730 std::array<MInt, nDim> endBuffer{0};
4731 std::array<MInt, nDim> sizeBuffer{0};
4734 stepBuffer[rcv->orderInfo[j]] = rcv->stepInfo[j];
4738 endBuffer[j] = size[j] - 1;
4739 sizeBuffer[rcv->orderInfo[j]] = size[j];
4740 if(stepBuffer[j] < 0) {
4741 std::swap(startBuffer[j], endBuffer[j]);
4746 auto* orderInfo = rcv->orderInfo.data();
4747 auto* startInfoCells = rcv->startInfoCells.data();
4748 auto* cellBuffer = rcv->cellBuffer.get();
4749 auto* variables = &(rcv->variables[0]);
4750 for(
MInt var = 0; var < rcv->noVars; var++) {
4751 maia::parallelFor<true, nDim>(begin, end, [=](
const MInt& i,
const MInt& j,
const MInt& k) {
4752 std::array<MInt, nDim> start{};
4753 start[orderInfo[0]] = startBuffer[0] + (i - startInfoCells[0]) * stepBuffer[0];
4754 start[orderInfo[1]] = startBuffer[1] + (j - startInfoCells[1]) * stepBuffer[1];
4755 start[orderInfo[2]] = startBuffer[2] + (k - startInfoCells[2]) * stepBuffer[2];
4757 const MInt bufferId = var * totalSize + start[0] + (start[1] + start[2] * sizeBuffer[1]) * sizeBuffer[0];
4759 variables[var][cellId] = cellBuffer[bufferId];
4766 std::vector<MPI_Request> sndRequests;
4767 std::vector<MPI_Request> rcvRequests;
4768 std::vector<MPI_Status> sndStatus;
4769 std::vector<MPI_Status> rcvStatus;
4777 sndStatus.resize(sndRequests.size());
4778 MPI_Waitall(sndRequests.size(), &sndRequests[0], &sndStatus[0], AT_);
4779 rcvStatus.resize(rcvRequests.size());
4780 MPI_Waitall(rcvRequests.size(), &rcvRequests[0], &rcvStatus[0], AT_);
4790 for(
MInt var = 0; var <
PV->noVariables; var++) {
4791 for(
MInt k = snd->startInfoCells[2]; k < snd->endInfoCells[2]; k++) {
4792 for(
MInt j = snd->startInfoCells[1]; j < snd->endInfoCells[1]; j++) {
4793 for(
MInt i = snd->startInfoCells[0]; i < snd->endInfoCells[0]; i++) {
4803 for(
MInt var = 0; var <
nDim; var++) {
4804 for(
MInt k = snd->startInfoCells[2]; k < snd->endInfoCells[2]; k++) {
4805 for(
MInt j = snd->startInfoCells[1]; j < snd->endInfoCells[1]; j++) {
4806 for(
MInt i = snd->startInfoCells[0]; i < snd->endInfoCells[0]; i++) {
4808 snd->cellBuffer[pos] =
m_cells->
fq[
FQ->VORTICITY[var]][cellId];
4820 MPI_Request request{};
4822 const MInt err =
MPI_Isend((
void*)&snd->cellBuffer[0], snd->cellBufferSize, MPI_DOUBLE, snd->nghbrId, tag,
4824 sndRequests.push_back(request);
4825 if(err) cout <<
"rank " <<
domainId() <<
" sending throws error " << endl;
4831 MPI_Request request{};
4832 const MInt tag = rcv->nghbrId + (rcv->tagHelper) *
noDomains();
4833 const MInt err =
MPI_Irecv((
void*)&rcv->cellBuffer[0], rcv->cellBufferSize, MPI_DOUBLE, rcv->nghbrId, tag,
4835 rcvRequests.push_back(request);
4836 if(err) cout <<
"rank " <<
domainId() <<
" sending throws error " << endl;
4844 for(
MInt var = 0; var <
PV->noVariables; var++) {
4845 for(
MInt k = rcv->startInfoCells[2]; k < rcv->endInfoCells[2]; k++) {
4846 for(
MInt j = rcv->startInfoCells[1]; j < rcv->endInfoCells[1]; j++) {
4847 for(
MInt i = rcv->startInfoCells[0]; i < rcv->endInfoCells[0]; i++) {
4857 for(
MInt var = 0; var < (2 *
nDim - 3); var++) {
4858 for(
MInt k = rcv->startInfoCells[2]; k < rcv->endInfoCells[2]; k++) {
4859 for(
MInt j = rcv->startInfoCells[1]; j < rcv->endInfoCells[1]; j++) {
4860 for(
MInt i = rcv->startInfoCells[0]; i < rcv->endInfoCells[0]; i++) {
4875 MInt allCellsK =
m_grid->getBlockNoCells(0, 0);
4885 MInt noVars =
PV->noVariables;
4887 noVars += (2 *
nDim - 3);
4901 std::vector<std::unique_ptr<StructuredComm<nDim>>> gcSndComm;
4902 std::vector<std::unique_ptr<StructuredComm<nDim>>> gcRcvComm;
4904 MFloat*
const*
const varPtr = &variables[0];
4906 m_windowInfo->createCommunicationExchangeFlags(gcSndComm, gcRcvComm, (
MInt)variables.size(), varPtr);
4913 MInt cellId, cellIdAdj1, cellIdAdj2;
4914 const MInt noVars = variables.size();
4922 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4923 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4930 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4931 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4945 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4946 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4953 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4954 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4968 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4969 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4976 for(
MInt varPos = 0; varPos < noVars; varPos++) {
4977 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4992 for(
MInt dim = 0; dim < 3; dim++) {
5019template <MBool twoEqRans>
5023 constexpr MFloat rPrTurb = F1 / 0.9;
5026 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
5047 T[I] =
m_gamma * p[I] / rho[I];
5048 muLam[I] = SUTHERLANDLAW(T[I]);
5062 const MFloat cornerMetrics[9] = {
5068 const MFloat dudxi = F1B4 * (u[IPJPKP] + u[IPJPK] + u[IPJKP] + u[IPJK] - u[IJPKP] - u[IJPK] - u[IJKP] - u[IJK]);
5069 const MFloat dudet = F1B4 * (u[IPJPKP] + u[IJPKP] + u[IPJPK] + u[IJPK] - u[IPJKP] - u[IJKP] - u[IPJK] - u[IJK]);
5070 const MFloat dudze = F1B4 * (u[IPJPKP] + u[IJPKP] + u[IPJKP] + u[IJKP] - u[IPJPK] - u[IJPK] - u[IPJK] - u[IJK]);
5072 const MFloat dvdxi = F1B4 * (v[IPJPKP] + v[IPJPK] + v[IPJKP] + v[IPJK] - v[IJPKP] - v[IJPK] - v[IJKP] - v[IJK]);
5073 const MFloat dvdet = F1B4 * (v[IPJPKP] + v[IJPKP] + v[IPJPK] + v[IJPK] - v[IPJKP] - v[IJKP] - v[IPJK] - v[IJK]);
5074 const MFloat dvdze = F1B4 * (v[IPJPKP] + v[IJPKP] + v[IPJKP] + v[IJKP] - v[IPJPK] - v[IJPK] - v[IPJK] - v[IJK]);
5076 const MFloat dwdxi = F1B4 * (w[IPJPKP] + w[IPJPK] + w[IPJKP] + w[IPJK] - w[IJPKP] - w[IJPK] - w[IJKP] - w[IJK]);
5077 const MFloat dwdet = F1B4 * (w[IPJPKP] + w[IJPKP] + w[IPJPK] + w[IJPK] - w[IPJKP] - w[IJKP] - w[IPJK] - w[IJK]);
5078 const MFloat dwdze = F1B4 * (w[IPJPKP] + w[IJPKP] + w[IPJKP] + w[IJKP] - w[IPJPK] - w[IJPK] - w[IPJK] - w[IJK]);
5080 const MFloat dTdxi = F1B4 * (T[IPJPKP] + T[IPJPK] + T[IPJKP] + T[IPJK] - T[IJPKP] - T[IJPK] - T[IJKP] - T[IJK]);
5081 const MFloat dTdet = F1B4 * (T[IPJPKP] + T[IJPKP] + T[IPJPK] + T[IJPK] - T[IPJKP] - T[IJKP] - T[IPJK] - T[IJK]);
5082 const MFloat dTdze = F1B4 * (T[IPJPKP] + T[IJPKP] + T[IPJKP] + T[IJKP] - T[IPJPK] - T[IJPK] - T[IPJK] - T[IJK]);
5084 const MFloat uAvg = F1B8 * (u[IPJPKP] + u[IJPKP] + u[IJPK] + u[IPJPK] + u[IPJKP] + u[IJKP] + u[IJK] + u[IPJK]);
5085 const MFloat vAvg = F1B8 * (v[IPJPKP] + v[IJPKP] + v[IJPK] + v[IPJPK] + v[IPJKP] + v[IJKP] + v[IJK] + v[IPJK]);
5086 const MFloat wAvg = F1B8 * (w[IPJPKP] + w[IJPKP] + w[IJPK] + w[IPJPK] + w[IPJKP] + w[IJKP] + w[IJK] + w[IPJK]);
5088 const MFloat muLamAvg = F1B8
5089 * (muLam[IPJPKP] + muLam[IJPKP] + muLam[IJPK] + muLam[IPJPK] + muLam[IPJKP] + muLam[IJKP]
5090 + muLam[IJK] + muLam[IPJK]);
5094 const MFloat muTurbAvg = F1B8
5095 * (muTurb[IPJPKP] + muTurb[IJPKP] + muTurb[IJPK] + muTurb[IPJPK] + muTurb[IPJKP]
5096 + muTurb[IJKP] + muTurb[IJK] + muTurb[IPJK]);
5102 * (rho[IPJPKP] * TKE[IPJPKP] + rho[IJPKP] * TKE[IJPKP] + rho[IJPK] * TKE[IJPK] + rho[IPJPK] * TKE[IPJPK]
5103 + rho[IPJKP] * TKE[IPJKP] + rho[IJKP] * TKE[IJKP] + rho[IJK] * TKE[IJK] + rho[IPJK] * TKE[IPJK]);
5111 * (dudxi * cornerMetrics[
xsd * 3 +
xsd] + dudet * cornerMetrics[
ysd * 3 +
xsd]
5112 + dudze * cornerMetrics[
zsd * 3 +
xsd])
5116 * (dvdxi * cornerMetrics[
xsd * 3 +
ysd] + dvdet * cornerMetrics[
ysd * 3 +
ysd]
5117 + dvdze * cornerMetrics[
zsd * 3 +
ysd])
5121 * (dwdxi * cornerMetrics[
xsd * 3 +
zsd] + dwdet * cornerMetrics[
ysd * 3 +
zsd]
5122 + dwdze * cornerMetrics[
zsd * 3 +
zsd]);
5128 MFloat tau2 = dudxi * cornerMetrics[
xsd * 3 +
ysd] + dudet * cornerMetrics[
ysd * 3 +
ysd]
5129 + dudze * cornerMetrics[
zsd * 3 +
ysd] +
5131 dvdxi * cornerMetrics[
xsd * 3 +
xsd] + dvdet * cornerMetrics[
ysd * 3 +
xsd]
5132 + dvdze * cornerMetrics[
zsd * 3 +
xsd];
5138 MFloat tau3 = dudxi * cornerMetrics[
xsd * 3 +
zsd] + dudet * cornerMetrics[
ysd * 3 +
zsd]
5139 + dudze * cornerMetrics[
zsd * 3 +
zsd] +
5141 dwdxi * cornerMetrics[
xsd * 3 +
xsd] + dwdet * cornerMetrics[
ysd * 3 +
xsd]
5142 + dwdze * cornerMetrics[
zsd * 3 +
xsd];
5150 * (dvdxi * cornerMetrics[
xsd * 3 +
ysd] + dvdet * cornerMetrics[
ysd * 3 +
ysd]
5151 + dvdze * cornerMetrics[
zsd * 3 +
ysd])
5155 * (dudxi * cornerMetrics[
xsd * 3 +
xsd] + dudet * cornerMetrics[
ysd * 3 +
xsd]
5156 + dudze * cornerMetrics[
zsd * 3 +
xsd])
5160 * (dwdxi * cornerMetrics[
xsd * 3 +
zsd] + dwdet * cornerMetrics[
ysd * 3 +
zsd]
5161 + dwdze * cornerMetrics[
zsd * 3 +
zsd]);
5167 MFloat tau5 = dvdxi * cornerMetrics[
xsd * 3 +
zsd] + dvdet * cornerMetrics[
ysd * 3 +
zsd]
5168 + dvdze * cornerMetrics[
zsd * 3 +
zsd] +
5170 dwdxi * cornerMetrics[
xsd * 3 +
ysd] + dwdet * cornerMetrics[
ysd * 3 +
ysd]
5171 + dwdze * cornerMetrics[
zsd * 3 +
ysd];
5179 * (dwdxi * cornerMetrics[
xsd * 3 +
zsd] + dwdet * cornerMetrics[
ysd * 3 +
zsd]
5180 + dwdze * cornerMetrics[
zsd * 3 +
zsd])
5184 * (dudxi * cornerMetrics[
xsd * 3 +
xsd] + dudet * cornerMetrics[
ysd * 3 +
xsd]
5185 + dudze * cornerMetrics[
zsd * 3 +
xsd])
5189 * (dvdxi * cornerMetrics[
xsd * 3 +
ysd] + dvdet * cornerMetrics[
ysd * 3 +
ysd]
5190 + dvdze * cornerMetrics[
zsd * 3 +
ysd]);
5193 const MFloat dTdx = dTdxi * cornerMetrics[
xsd * 3 +
xsd] + dTdet * cornerMetrics[
ysd * 3 +
xsd]
5194 + dTdze * cornerMetrics[
zsd * 3 +
xsd];
5196 const MFloat dTdy = dTdxi * cornerMetrics[
xsd * 3 +
ysd] + dTdet * cornerMetrics[
ysd * 3 +
ysd]
5197 + dTdze * cornerMetrics[
zsd * 3 +
ysd];
5199 const MFloat dTdz = dTdxi * cornerMetrics[
xsd * 3 +
zsd] + dTdet * cornerMetrics[
ysd * 3 +
zsd]
5200 + dTdze * cornerMetrics[
zsd * 3 +
zsd];
5203 const MFloat mueOverRe = rRe * fJac * (muLamAvg + muTurbAvg);
5204 tau1 = mueOverRe * tau1 + TKEcorner;
5207 tau4 = mueOverRe * tau4 + TKEcorner;
5209 tau6 = mueOverRe * tau6 + TKEcorner;
5211 const MFloat muCombined = FgammaMinusOne * rRe * fJac * (rPrLam * muLamAvg + rPrTurb * muTurbAvg);
5213 const MFloat qx = muCombined * dTdx + uAvg * tau1 + vAvg * tau2 + wAvg * tau3;
5214 const MFloat qy = muCombined * dTdy + uAvg * tau2 + vAvg * tau4 + wAvg * tau5;
5215 const MFloat qz = muCombined * dTdz + uAvg * tau3 + vAvg * tau5 + wAvg * tau6;
5219 tau1 * cornerMetrics[
xsd * 3 +
xsd] + tau2 * cornerMetrics[
xsd * 3 +
ysd] + tau3 * cornerMetrics[
xsd * 3 +
zsd];
5222 tau2 * cornerMetrics[
xsd * 3 +
xsd] + tau4 * cornerMetrics[
xsd * 3 +
ysd] + tau5 * cornerMetrics[
xsd * 3 +
zsd];
5225 tau3 * cornerMetrics[
xsd * 3 +
xsd] + tau5 * cornerMetrics[
xsd * 3 +
ysd] + tau6 * cornerMetrics[
xsd * 3 +
zsd];
5228 qx * cornerMetrics[
xsd * 3 +
xsd] + qy * cornerMetrics[
xsd * 3 +
ysd] + qz * cornerMetrics[
xsd * 3 +
zsd];
5232 tau1 * cornerMetrics[
ysd * 3 +
xsd] + tau2 * cornerMetrics[
ysd * 3 +
ysd] + tau3 * cornerMetrics[
ysd * 3 +
zsd];
5235 tau2 * cornerMetrics[
ysd * 3 +
xsd] + tau4 * cornerMetrics[
ysd * 3 +
ysd] + tau5 * cornerMetrics[
ysd * 3 +
zsd];
5238 tau3 * cornerMetrics[
ysd * 3 +
xsd] + tau5 * cornerMetrics[
ysd * 3 +
ysd] + tau6 * cornerMetrics[
ysd * 3 +
zsd];
5241 qx * cornerMetrics[
ysd * 3 +
xsd] + qy * cornerMetrics[
ysd * 3 +
ysd] + qz * cornerMetrics[
ysd * 3 +
zsd];
5245 tau1 * cornerMetrics[
zsd * 3 +
xsd] + tau2 * cornerMetrics[
zsd * 3 +
ysd] + tau3 * cornerMetrics[
zsd * 3 +
zsd];
5248 tau2 * cornerMetrics[
zsd * 3 +
xsd] + tau4 * cornerMetrics[
zsd * 3 +
ysd] + tau5 * cornerMetrics[
zsd * 3 +
zsd];
5251 tau3 * cornerMetrics[
zsd * 3 +
xsd] + tau5 * cornerMetrics[
zsd * 3 +
ysd] + tau6 * cornerMetrics[
zsd * 3 +
zsd];
5254 qx * cornerMetrics[
zsd * 3 +
xsd] + qy * cornerMetrics[
zsd * 3 +
ysd] + qz * cornerMetrics[
zsd * 3 +
zsd];
5264 for(
MInt var = 0; var < 4; var++) {
5266 maia::parallelFor<true, nDim>(beginIM1,
endM2(), [=](
const MInt& i,
const MInt& j,
const MInt& k) {
5272 vflux[0][IJK] = F1B4 * (eflux[var][IJK] + eflux[var][IJKM] + eflux[var][IJMK] + eflux[var][IJMKM]);
5277 maia::parallelFor<true, nDim>(beginJM1,
endM2(), [=](
const MInt& i,
const MInt& j,
const MInt& k) {
5283 vflux[1][IJK] = F1B4 * (fflux[var][IJK] + fflux[var][IJKM] + fflux[var][IMJK] + fflux[var][IMJKM]);
5287 maia::parallelFor<true, nDim>(beginKM1,
endM2(), [=](
const MInt& i,
const MInt& j,
const MInt& k) {
5293 vflux[2][IJK] = F1B4 * (gflux[var][IJK] + gflux[var][IMJK] + gflux[var][IJMK] + gflux[var][IMJMK]);
5302 vflux[0][IJK] - vflux[0][IMJK] + vflux[1][IJK] - vflux[1][IJMK] + vflux[2][IJK] - vflux[2][IJKM];
5312 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
5325 MInt start[3], end[3], nghbr[20];
5335 if(len1[j] != 0) totalCells *= len1[j];
5338 for(
MInt n = 0; n < 3; ++n) {
5350 MFloat u[20], v[20], w[20], T[20];
5351 MFloat U, V, W, t, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, dTdx, dTdy, dTdz;
5353 MFloat tau1, tau2, tau3, tau4, tau5, tau6;
5354 MFloat mueOverRe, mue, mueH;
5355 for(
MInt kk = start[2]; kk < end[2]; ++kk) {
5356 for(
MInt jj = start[1]; jj < end[1]; ++jj) {
5357 for(
MInt ii = start[0]; ii < end[0]; ++ii) {
5359 MInt temp[3] = {0, 0, 0};
5363 const MFloat cornerMetrics[9] = {
5370 nghbr[count++] =
cellIndex(ii + temp[0], jj + temp[1], kk + temp[2]);
5374 nghbr[count++] =
cellIndex(ii + change[0], jj + change[1], kk + change[2]);
5375 nghbr[count++] =
cellIndex(ii + temp[0] + change[0], jj + temp[1] + change[1], kk + temp[2] + change[2]);
5379 cout <<
"what the hell! it is wrong!!!" << endl;
5383 u[m] = rhou[nghbr[m]] / rho[nghbr[m]];
5384 v[m] = rhov[nghbr[m]] / rho[nghbr[m]];
5385 w[m] = rhow[nghbr[m]] / rho[nghbr[m]];
5386 T[m] = (
m_gamma * gammaMinusOne
5387 * (rhoE[nghbr[m]] - F1B2 * rho[nghbr[m]] * (
POW2(u[m]) +
POW2(v[m]) +
POW2(w[m]))))
5408 MInt id2 = ii - start[0] + ((jj - start[1]) + (kk - start[2]) * len1[1]) * len1[0];
5410 for(
MInt n = 0; n < count; n++) {
5411 MInt ID = id2 * count + n;
5433 tau1 = 2 * dudx - 2 / 3 * (dudx + dvdy + dwdz);
5436 tau4 = 2 * dvdy - 2 / 3 * (dudx + dvdy + dwdz);
5438 tau6 = 2 * dwdz - 2 / 3 * (dudx + dvdy + dwdz);
5440 mue = SUTHERLANDLAW(t);
5441 mueOverRe = mue * rRe;
5448 mueH = FgammaMinusOne * mueOverRe * rPr;
5450 qx = mueH * dTdx + U * tau1 + V * tau2 + W * tau3;
5451 qy = mueH * dTdy + U * tau2 + V * tau4 + W * tau5;
5452 qz = mueH * dTdz + U * tau3 + V * tau5 + W * tau6;
5456 eflux[0][IJK] = tau1 * cornerMetrics[
xsd * 3 +
xsd] + tau2 * cornerMetrics[
xsd * 3 +
ysd]
5457 + tau3 * cornerMetrics[
xsd * 3 +
zsd];
5459 eflux[1][IJK] = tau2 * cornerMetrics[
xsd * 3 +
xsd] + tau4 * cornerMetrics[
xsd * 3 +
ysd]
5460 + tau5 * cornerMetrics[
xsd * 3 +
zsd];
5462 eflux[2][IJK] = tau3 * cornerMetrics[
xsd * 3 +
xsd] + tau5 * cornerMetrics[
xsd * 3 +
ysd]
5463 + tau6 * cornerMetrics[
xsd * 3 +
zsd];
5465 eflux[3][IJK] = qx * cornerMetrics[
xsd * 3 +
xsd] + qy * cornerMetrics[
xsd * 3 +
ysd]
5466 + qz * cornerMetrics[
xsd * 3 +
zsd];
5469 fflux[0][IJK] = tau1 * cornerMetrics[
ysd * 3 +
xsd] + tau2 * cornerMetrics[
ysd * 3 +
ysd]
5470 + tau3 * cornerMetrics[
ysd * 3 +
zsd];
5472 fflux[1][IJK] = tau2 * cornerMetrics[
ysd * 3 +
xsd] + tau4 * cornerMetrics[
ysd * 3 +
ysd]
5473 + tau5 * cornerMetrics[
ysd * 3 +
zsd];
5475 fflux[2][IJK] = tau3 * cornerMetrics[
ysd * 3 +
xsd] + tau5 * cornerMetrics[
ysd * 3 +
ysd]
5476 + tau6 * cornerMetrics[
ysd * 3 +
zsd];
5478 fflux[3][IJK] = qx * cornerMetrics[
ysd * 3 +
xsd] + qy * cornerMetrics[
ysd * 3 +
ysd]
5479 + qz * cornerMetrics[
ysd * 3 +
zsd];
5482 gflux[0][IJK] = tau1 * cornerMetrics[
zsd * 3 +
xsd] + tau2 * cornerMetrics[
zsd * 3 +
ysd]
5483 + tau3 * cornerMetrics[
zsd * 3 +
zsd];
5485 gflux[1][IJK] = tau2 * cornerMetrics[
zsd * 3 +
xsd] + tau4 * cornerMetrics[
zsd * 3 +
ysd]
5486 + tau5 * cornerMetrics[
zsd * 3 +
zsd];
5488 gflux[2][IJK] = tau3 * cornerMetrics[
zsd * 3 +
xsd] + tau5 * cornerMetrics[
zsd * 3 +
ysd]
5489 + tau6 * cornerMetrics[
zsd * 3 +
zsd];
5491 gflux[3][IJK] = qx * cornerMetrics[
zsd * 3 +
xsd] + qy * cornerMetrics[
zsd * 3 +
ysd]
5492 + qz * cornerMetrics[
zsd * 3 +
zsd];
5513 mTerm(1,
"Porous stuff for 3D RANS not implemented yet!");
5520 for(
MInt dim = 0; dim <
nDim; ++dim) {
5521 velAbs +=
POW2(pvars[
PV->VV[dim]][IJK]);
5523 velAbs = sqrt(velAbs);
5524 const MFloat rDa = 1.0 / Da[IJK];
5526 for(
MInt dim = 0; dim <
nDim; ++dim) {
5528 -(rRe0 * rDa * por[IJK] * muLam[IJK] + porPOW2 * sqrt(rDa) * cf[IJK] * velAbs * pvars[
PV->RHO][IJK])
5558 MFloat epsilon = pow(10.0, -10.0);
5562 MFloat maxResidual1 = F0;
5563 MInt maxResIndex[3];
5565 MFloat maxResidualOrg = F0;
5566 MFloat localMaxResidual = F0;
5567 MFloat localAvrgResidual = F0;
5568 MFloat accumAvrgResidual = F0;
5569 MFloat globalMaxResidual = F0;
5572 for(
MInt dim = 0; dim <
nDim; dim++) {
5573 maxResIndex[dim] = F0;
5584 if(tmpResidual > maxResidual1) {
5588 maxResidual1 = tmpResidual;
5595 localMaxResidual = maxResidual1;
5600 MFloat localTotalEnergy = F0;
5601 MFloat globalTotalEnergy = F0;
5602 MFloat globalPressure = F0;
5603 MFloat localPressure = F0;
5607 "localTotalEnergy",
"globalTotalEnergy");
5608 globalTotalEnergy /= ((16.0 * pow(4 * atan(1), 3.0)));
5612 globalPressure /= ((16.0 * pow(4 * atan(1), 3.0)));
5620 "localAvrgResidual",
"accumAvrgResidual");
5622 "localMaxResidual",
"globalMaxResidual");
5625 maxResidualOrg = globalMaxResidual;
5641 if(
approx(localMaxResidual, maxResidualOrg,
5646 f_residual = fopen(
"./Residual",
"a+");
5649 fprintf(f_residual,
"#iter, physTime, time, dT, wLoad, avrgRes, maxRes, blockId, i, j, k ");
5656 f_residual = fopen(
"./Residual",
"a+");
5660 "#iter, physTime, time, dT, wLoad, avrgRes, maxRes, blockId, i, j, k, k_mean, p_mean, pProbe ");
5673 cerr <<
"Solution diverged, average residual is nan " << endl;
5674 m_log <<
"Solution diverged, average residual is nan " << endl;
5676 mTerm(1, AT_,
"Solution diverged, average residual is nan ");
5690 if(
approx(localMaxResidual, maxResidualOrg,
5694 f_residual = fopen(
"./Residual",
"a+");
5697 fprintf(f_residual,
" %f",
m_time);
5701 fprintf(f_residual,
" %1.10e", globalMaxResidual);
5706 fprintf(f_residual,
"\n");
5722 f_residual = fopen(
"./Residual",
"a+");
5725 fprintf(f_residual,
" %f",
m_time);
5729 fprintf(f_residual,
" %1.10e", globalMaxResidual);
5734 fprintf(f_residual,
" %1.10e", globalTotalEnergy);
5735 fprintf(f_residual,
" %1.10e", globalPressure);
5736 fprintf(f_residual,
" %1.10e", dissip);
5737 fprintf(f_residual,
"\n");
5767 MFloat localPressure = F0;
5776 return localPressure;
5800 const MFloat pi = 4.0 * atan(1);
5808 const MFloat y_max = 10.0;
5810 const MFloat amp =
PV->UInfinity / (F2 * pi * frequency);
5816 const MFloat x =
m_grid->m_initCoordinates[0][pointId];
5818 m_grid->m_coordinates[0][pointId] =
5819 x + (1 -
m_grid->m_coordinates[1][pointId] / y_max) * amp * sin(2 * pi * t * frequency);
5820 m_grid->m_velocity[0][pointId] = (1 -
m_grid->m_coordinates[1][pointId] / y_max) * amp * 2 * pi * frequency
5821 * cos(2 * pi * t * frequency);
5830 const MFloat beta = 4.14l;
5836 const MFloat x2 = y_max * (ver - 11.75l);
5837 const MFloat x3 = y_max * (ver - 9.25l);
5838 const MFloat x4 = y_max * (ver - 1.25);
5839 const MFloat x5 = y_max * (ver + 1.25);
5840 const MFloat omega =
PV->UInfinity * StrNum * F2 * pi / y_max;
5841 const MFloat h = F1B2 * 0.38l * (F1 - cos(omega * t));
5842 const MFloat hvel = F1B2 * 0.38l * omega * sin(omega * t);
5852 if(
m_grid->m_coordinates[0][pointId] > x2 &&
m_grid->m_coordinates[0][pointId] < x3) {
5853 g = (1.0 + tanhl(beta * (
m_grid->m_coordinates[0][pointId] - (x2 + x3) / 2.0l) / y_max)) / 2.0;
5854 }
else if(
m_grid->m_coordinates[0][pointId] >= x3 &&
m_grid->m_coordinates[0][pointId] < x4) {
5856 }
else if(
m_grid->m_coordinates[0][pointId] >= x4 &&
m_grid->m_coordinates[0][pointId] < x5) {
5857 g = (1.0 - tanhl(beta * (
m_grid->m_coordinates[0][pointId] - (x4 + x5) / 2.0l) / y_max)) / 2.0;
5860 m_grid->m_coordinates[1][pointId] =
y * (F1 - h * g);
5861 m_grid->m_velocity[1][pointId] = -
y * hvel * g;
5874 const MFloat x =
m_grid->m_initCoordinates[0][pointId];
5878 m_grid->m_acceleration[0][pointId] = F0;
5887 const MFloat beta = 16.0l;
5890 const MFloat y_max = 1.0l;
5893 const MFloat x2 = y_max * (ver + 0.1l);
5894 const MFloat x3 = y_max * (ver + 0.5l);
5895 const MFloat x4 = y_max * (ver + 0.5l);
5896 const MFloat x5 = y_max * (ver + 0.9l);
5898 const MFloat omega =
m_Ma * sqrt(
PV->TInfinity) * StrNum * 2.0l * pi / y_max;
5899 const MFloat h = F1B2 * 0.35l * (F1 - cos(omega * t));
5900 const MFloat hvel = F1B2 * 0.35l * (omega * sin(omega * t));
5907 const MFloat x =
m_grid->m_initCoordinates[0][pointId];
5911 if(
y > x2 && y < x5 && x > x2 && x < x5) {
5912 g = ((
y < x3) ? (1.0l + tanhl(beta * (
y - (x2 + x3) / 2.0l) / y_max)) / 2.0l
5913 : ((
y < x4) ? 1.0l : (1.0l - tanhl(beta * (
y - (x4 + x5) / 2.0l) / y_max)) / 2.0l));
5915 m_grid->m_coordinates[0][pointId] = x * (F1 - h * g * (F1 - x));
5916 m_grid->m_velocity[0][pointId] = x * (-hvel * g * (F1 - x));
5919 if(x > x2 && x < x5 && y > x2 &&
y < x5) {
5920 g = ((x < x3) ? (1.0l + tanhl(beta * (x - (x2 + x3) / 2.0l) / y_max)) / 2.0l
5921 : ((x < x4) ? 1.0l : (1.0l - tanhl(beta * (x - (x4 + x5) / 2.0l) / y_max)) / 2.0l));
5923 m_grid->m_coordinates[1][pointId] =
y * (F1 - h * g * (F1 -
y));
5924 m_grid->m_velocity[1][pointId] =
y * (-hvel * g * (F1 -
y));
5935 const MFloat rotSin = sin(angle);
5944 MFloat fadeInFactor = 0.0;
5945 MFloat fadeInFactorPrime = 0.0;
5946 MFloat fadeInFactorPrimePrime = 0.0;
5950 fadeInFactorPrimePrime =
5954 fadeInFactorPrime = 0.0;
5955 fadeInFactorPrimePrime = 0.0;
5967 const MFloat xInit =
m_grid->m_initCoordinates[0][pointId];
5968 const MFloat zInit =
m_grid->m_initCoordinates[2][pointId];
5969 const MFloat yInit =
m_grid->m_initCoordinates[1][pointId];
5971 MFloat transitionFactor = F0;
5973 transitionFactor = F0;
5977 transitionFactor = F1;
5981 transitionFactor = F0;
5984 MFloat yTransitionFactor = F1;
5986 yTransitionFactor = F1;
5990 yTransitionFactor = F0;
5995 const MFloat func = transitionFactor * yTransitionFactor
6000 const MFloat funcPrimePrime = -transitionFactor * yTransitionFactor
6004 m_grid->m_coordinates[1][pointId] = func * fadeInFactor + yInit;
6005 m_grid->m_velocity[1][pointId] = func * fadeInFactorPrime + funcPrime * fadeInFactor;
6006 m_grid->m_acceleration[1][pointId] =
6007 funcPrimePrime * fadeInFactor + 2 * funcPrime * fadeInFactorPrime + func * fadeInFactorPrimePrime;
6022 const MFloat timeRelaxation = 50.0;
6024 if(t_offset < timeRelaxation) {
6025 fadeInFactor = (1.0 - cos(t_offset / timeRelaxation * pi)) * F1B2;
6038 const MFloat yInit =
m_grid->m_initCoordinates[1][pointId];
6039 const MFloat zInit =
m_grid->m_initCoordinates[2][pointId];
6044 }
else if(yInit > F0 && yInit < F1) {
6045 yRelaxation = F1 - yInit;
6046 }
else if(yInit > F1 && yInit < F2) {
6047 yRelaxation = yInit - F1;
6053 m_grid->m_coordinates[1][pointId] =
6058 m_grid->m_velocity[1][pointId] = fadeInFactor
6062 m_grid->m_coordinates[1][pointId] =
6067 m_grid->m_velocity[1][pointId] = -fadeInFactor
6087 MFloat fadeInFactor = 0.0;
6088 MFloat fadeInFactorPrime = 0.0;
6089 MFloat fadeInFactorPrimePrime = 0.0;
6093 fadeInFactorPrimePrime =
6097 fadeInFactorPrime = 0.0;
6098 fadeInFactorPrimePrime = 0.0;
6110 const MFloat xInit =
m_grid->m_initCoordinates[0][pointId];
6111 const MFloat yInit =
m_grid->m_initCoordinates[1][pointId];
6113 MFloat transitionFactor = F0;
6115 transitionFactor = F0;
6119 transitionFactor = F1;
6123 transitionFactor = F0;
6126 MFloat yTransitionFactor = F1;
6128 yTransitionFactor = F1;
6132 yTransitionFactor = F0;
6135 const MFloat func = transitionFactor * yTransitionFactor
6140 const MFloat funcPrimePrime = -transitionFactor * yTransitionFactor
6144 m_grid->m_coordinates[1][pointId] = func * fadeInFactor + yInit;
6145 m_grid->m_velocity[1][pointId] = func * fadeInFactorPrime + funcPrime * fadeInFactor;
6146 m_grid->m_acceleration[1][pointId] =
6147 funcPrimePrime * fadeInFactor + 2 * funcPrime * fadeInFactorPrime + func * fadeInFactorPrimePrime;
6165 const MFloat timeRelaxation = 1.0;
6167 if(t_offset < timeRelaxation) {
6168 fadeInFactor = (1.0 - cos(t_offset / timeRelaxation * pi)) * F1B2;
6177 const MInt myBlockId =
m_grid->getMyBlockId();
6179 if(myBlockId == 0) {
6181 MFloat transitionFactor = F0;
6183 const MInt globalPointId = offset + i;
6192 transitionFactor = F0;
6196 transitionFactor = F1;
6200 transitionFactor = F0;
6208 if(wallNormalVec[1] < 0.0) {
6217 const MFloat xInit =
m_grid->m_initCoordinates[0][pIJK];
6218 const MFloat yInit =
m_grid->m_initCoordinates[1][pIJK];
6219 const MFloat zInit =
m_grid->m_initCoordinates[2][pIJK];
6222 MFloat yTransitionFactor = F1;
6223 const MFloat normalDist = sqrt(
POW2(xInit - xWall) +
POW2(yInit - yWall));
6225 yTransitionFactor = F1;
6229 yTransitionFactor = F0;
6232 const MFloat normalDisplacement =
6234 * (amp0 * yTransitionFactor * transitionFactor
6241 m_grid->m_coordinates[0][pIJK] = xInit + normalDisplacement * wallNormalVec[0];
6242 m_grid->m_coordinates[1][pIJK] = yInit + normalDisplacement * wallNormalVec[1];
6244 m_grid->m_velocity[0][pIJK] = normalVel * wallNormalVec[0];
6245 m_grid->m_velocity[1][pIJK] = normalVel * wallNormalVec[1];
6260 const MFloat timeRelaxation = 1.0;
6262 if(t_offset < timeRelaxation) {
6263 fadeInFactor = (1.0 - cos(t_offset / timeRelaxation * pi)) * F1B2;
6282 MFloat spaceTransition = F1;
6285 spaceTransition = F0;
6286 }
else if(r >= 1.0 && r <= 41.0) {
6287 spaceTransition = fabs(r - 1.0) / 40.0;
6289 spaceTransition = F1;
6292 m_grid->m_coordinates[1][pIJK] =
6295 m_grid->m_velocity[1][pIJK] =
6296 (1.0 - spaceTransition)
6298 m_grid->m_acceleration[1][pIJK] =
6299 (1.0 - spaceTransition)
6320 const MFloat timeRelaxation = 1.0;
6322 if(t_offset < timeRelaxation) {
6323 fadeInFactor = (1.0 - cos(t_offset / timeRelaxation * pi)) * F1B2;
6332 const MInt myBlockId =
m_grid->getMyBlockId();
6334 if(myBlockId == 0) {
6336 MFloat transitionFactor = F0;
6338 const MInt globalPointId = offset + i;
6349 transitionFactor = F0;
6353 transitionFactor = F1;
6357 transitionFactor = F0;
6361 if(wallNormalVec[1] < 0.0) {
6364 transitionFactor = F0;
6369 transitionFactor = F1;
6374 transitionFactor = F0;
6381 const MFloat xInit =
m_grid->m_initCoordinates[0][pIJK];
6382 const MFloat yInit =
m_grid->m_initCoordinates[1][pIJK];
6384 MFloat yTransitionFactor = F1;
6385 const MFloat normalDist = sqrt(
POW2(xInit - xWall) +
POW2(yInit - yWall));
6387 yTransitionFactor = F1;
6391 yTransitionFactor = F0;
6394 const MFloat normalDisplacement =
6396 * (amp0 * yTransitionFactor * transitionFactor
6403 m_grid->m_coordinates[0][pIJK] = xInit + normalDisplacement * wallNormalVec[0];
6404 m_grid->m_coordinates[1][pIJK] = yInit + normalDisplacement * wallNormalVec[1];
6406 m_grid->m_velocity[0][pIJK] = normalVel * wallNormalVec[0];
6407 m_grid->m_velocity[1][pIJK] = normalVel * wallNormalVec[1];
6415 mTerm(1, AT_,
"Grid Moving Method not implemented!");
6423 m_grid->extrapolateGhostPointCoordinates();
6436 m_grid->computeCellCenterCoordinates();
6440 m_grid->computeMetrics();
6444 m_grid->computeJacobian();
6451 m_log <<
"Initializing moving grid methods..." << endl;
6459 m_grid->m_initCoordinates[isd][pointId] =
m_grid->m_coordinates[isd][pointId];
6508 mTerm(1, AT_,
"Property waveLengthPlus not specified in property file");
6528 mTerm(1, AT_,
"Property waveAmplitudePlus not specified in property file");
6547 mTerm(1, AT_,
"Property waveTimePlus not specified in property file");
6693 deltaS = pow(
m_Re, -7.0 / 8.0) * sqrt(2.0 / 0.024);
6694 cf = 0.024 * pow(
m_Re, -F1B4);
6697 deltaS = ((log(
m_Re)) / 0.384 + 4.127) /
m_Re;
6698 cf = 2.0 * pow((log(
m_Re) / 0.384 + 4.127), -2.0);
6700 const MFloat uTau = sqrt(
POW2(
PV->UInfinity) *
CV->rhoInfinity * cf * F1B2);
6715 m_log <<
"/////////////////// TRAVELING WAVE /////////////////////////////" << endl;
6716 m_log <<
"Re: " <<
m_Re <<
" c_f: " << cf <<
" u_tau: " << uTau <<
" deltaZ: " << deltaZ
6717 <<
" viscousUnit delta_v: " << deltaS << endl;
6722 m_log <<
"Max up/down speed: " << speedAmplitude << endl;
6724 m_log <<
"////////////////////////////////////////////////////////////////" << endl;
6762 mTerm(1, AT_,
"Property waveLength not specified in property file");
6781 mTerm(1, AT_,
"Property waveAmplitude not specified in property file");
6800 mTerm(1, AT_,
"Property waveTime not specified in property file");
6948 m_log <<
"/////////////////// TRAVELING WAVE /////////////////////////////" << endl;
6951 m_log <<
"Max up/down speed: " << speedAmplitude << endl;
6953 m_log <<
"////////////////////////////////////////////////////////////////" << endl;
6982 const MFloat mu8 = SUTHERLANDLAW(
PV->TInfinity);
7027 mTerm(1, AT_,
"Property waveLength not specified in property file");
7063 mTerm(1, AT_,
"Property waveAmplitude not specified in property file");
7082 mTerm(1, AT_,
"Property waveTime not specified in property file");
7208 m_log <<
"/////////////////// TRAVELING WAVE /////////////////////////////" << endl;
7211 m_log <<
"Speed amplitude: " << speedAmplitude << endl;
7213 m_log <<
"////////////////////////////////////////////////////////////////" << endl;
7403 const MInt myBlockId =
m_grid->getMyBlockId();
7404 const MInt blockId = 0;
7413 localCoords.fill(-99999.0);
7414 localNormalVec.fill(-99999.0);
7418 if(myBlockId == blockId) {
7423 const MInt localPointId = i;
7424 const MInt globalPointId = offset + i;
7430 MFloat normalVec[3] = {F0, F0, F0};
7431 for(
MInt dim = 0; dim <
nDim; dim++) {
7432 normalVec[dim] =
m_grid->m_coordinates[dim][pIJPK] -
m_grid->m_coordinates[dim][pIJK];
7436 const MFloat normalLength = sqrt(
POW2(normalVec[0]) +
POW2(normalVec[1]) +
POW2(normalVec[2]));
7437 for(
MInt dim = 0; dim <
nDim; dim++) {
7438 normalVec[dim] /= normalLength;
7441 localCoords(globalPointId) =
m_grid->m_coordinates[0][localPointId];
7445 localNormalVec(globalPointId) = normalVec[0];
7464 m_log <<
"/////////////////// TRAVELING WAVE /////////////////////////////" << endl;
7473 m_log <<
"Max up/down speed: " << speedAmplitude << endl;
7480 m_log <<
"////////////////////////////////////////////////////////////////" << endl;
7509 freqFactor = Context::getSolverProperty<MFloat>(
"oscFreqFactor",
m_solverId, AT_, &freqFactor);
7744 const MInt myBlockId =
m_grid->getMyBlockId();
7745 const MInt blockId = 0;
7754 localCoords.fill(-99999.0);
7755 localNormalVec.fill(-99999.0);
7759 if(myBlockId == blockId) {
7764 const MInt localPointId = i;
7765 const MInt globalPointId = offset + i;
7771 MFloat normalVec[3] = {F0, F0, F0};
7772 for(
MInt dim = 0; dim <
nDim; dim++) {
7773 normalVec[dim] =
m_grid->m_coordinates[dim][pIJPK] -
m_grid->m_coordinates[dim][pIJK];
7777 const MFloat normalLength = sqrt(
POW2(normalVec[0]) +
POW2(normalVec[1]) +
POW2(normalVec[2]));
7778 for(
MInt dim = 0; dim <
nDim; dim++) {
7779 normalVec[dim] /= normalLength;
7782 localCoords(globalPointId) =
m_grid->m_coordinates[0][localPointId];
7786 localNormalVec(globalPointId) = normalVec[0];
7805 m_log <<
"/////////////////// TRAVELING WAVE /////////////////////////////" << endl;
7812 m_log <<
"Max up/down speed: " << speedAmplitude << endl;
7823 m_log <<
"////////////////////////////////////////////////////////////////" << endl;
7827 stringstream errorMessage;
7828 errorMessage <<
"Grid moving method " <<
m_gridMovingMethod <<
" not implemented!" << endl;
7829 mTerm(1, AT_, errorMessage.str());
7860 m_log <<
"Initializing moving grid methods... DONE!" << endl;
7865 const MFloat pi = 4.0 * atan(1);
7878 MFloat fadeInFactor = 0.0;
7897 MFloat transitionFactor = F0;
7899 transitionFactor = F0;
7903 transitionFactor = F1;
7907 transitionFactor = F0;
7910 const MFloat force = fadeInFactor * transitionFactor
7933 MFloat fadeInFactor = 0.0;
7952 MFloat transitionFactor = F0;
7954 transitionFactor = F0;
7958 transitionFactor = F1;
7962 transitionFactor = F0;
7973 const MInt p1end = jj + (
m_grid->getMyBlockNoCells(0) - 1) *
m_grid->getMyBlockNoCells(1);
7978 amp = fstart + (pos - zzstart) / (zz - zzstart) * (f - fstart);
7980 const MInt p1 = jj + (
m_grid->getMyBlockNoCells(0) - 1) *
m_grid->getMyBlockNoCells(1);
7981 const MInt p1start = jj + 0 *
m_grid->getMyBlockNoCells(1);
7986 amp = f + (pos - zz) / (zzend - zz) * (fend - f);
7988 for(
MInt kk = 0; kk <
m_grid->getMyBlockNoCells(0) - 1; ++kk) {
7989 const MInt p1 = jj + kk *
m_grid->getMyBlockNoCells(1);
7990 const MInt p1p = jj + (kk + 1) *
m_grid->getMyBlockNoCells(1);
7995 if(zz <= pos && pos < zzp) {
7996 amp = f + (pos - zz) / (zzp - zz) * (fp - f);
8015 mTerm(1, AT_,
"Body Force Method not implemented!");
8022 m_log <<
"Initializing body force methods..." << endl;
8198 m_log <<
"/////////////////// TRAVELING WAVE BODY FORCE //////////////////" << endl;
8201 m_log <<
"Max up/down speed: " << speedAmplitude << endl;
8203 m_log <<
"////////////////////////////////////////////////////////////////" << endl;
8369 ParallelIo::size_type size[2] = {
m_grid->getMyBlockNoCells(0),
m_grid->getMyBlockNoCells(1)};
8377 stringstream pathStr;
8379 std::string path = pathStr.str();
8381 ParallelIo::size_type offset[2] = {0, 0};
8384 pio.readArray(&
m_waveForceY[0], path,
"y", 2, offset, size);
8385 pio.readArray(&
m_waveForceZ[0], path,
"z", 2, offset, size);
8399 m_log <<
"/////////////////// TRAVELING WAVE BODY FORCE //////////////////" << endl;
8402 m_log <<
"Max up/down speed: " << speedAmplitude << endl;
8404 m_log <<
"////////////////////////////////////////////////////////////////" << endl;
8408 "localDomainWidth",
"m_waveDomainWidth");
8414 stringstream errorMessage;
8415 errorMessage <<
"Grid moving method " <<
m_gridMovingMethod <<
" not implemented!" << endl;
8416 mTerm(1, AT_, errorMessage.str());
8420 m_log <<
"Initializing body force methods... DONE!" << endl;
8425 stringstream filename;
8462 for(
MInt dim = 0; dim <
nDim; ++dim) {
8487 for(
MInt dim = 0; dim <
nDim; ++dim) {
8519 const MInt localTotalId =
8523 for(
MInt dim = 0; dim <
nDim; dim++) {
8538 stringstream pathName;
8539 pathName <<
"box" <<
b;
8551 MInt hasCoordinates = 0;
8566 pio.
setAttribute(hasCoordinates,
"hasCoordinates", pathName.str());
8582 stringstream pathName;
8583 pathName <<
"box" <<
b;
8587 nDim, localOffset, localSize);
8599 ParallelIo::size_type localBoxPoints[3] = {0, 0, 0};
8600 ParallelIo::size_type localBoxOffset[3] = {0, 0, 0};
8608 pio.
writeArray(&empty, pathName.str(),
"x",
nDim, localBoxOffset, localBoxPoints);
8609 pio.
writeArray(&empty, pathName.str(),
"y",
nDim, localBoxOffset, localBoxPoints);
8610 pio.
writeArray(&empty, pathName.str(),
"z",
nDim, localBoxOffset, localBoxPoints);
8620 cout <<
"Loading BC2600 values..." << endl;
8624 MInt noCellsBC = bcCells[0] * bcCells[1] * bcCells[2];
8625 ParallelIo::size_type bcOffset[3] = {0, 0, 0};
8628 stringstream restartFileName;
8633 stringstream pathStr;
8634 pathStr <<
"/block" <<
m_blockId <<
"/bc2600";
8635 std::string path = pathStr.str();
8643 "tmpRestartVars(0");
8645 MInt startGC[3] = {0, 0, 0};
8646 MInt endGC[3] = {0, 0, 0};
8668 MInt cellIdBC = globalI + (globalJ + globalK * bcCells[1]) * bcCells[2];
8706 const MInt cellIdA2 =
8708 const MInt cellIdA1 =
8710 const MInt cellIdG1 =
8731 cout <<
"Loading BC2600 values... SUCCESSFUL!" << endl;
8742 cout <<
"Loading restart values 2601" << endl;
8744 ParallelIo::size_type bcCells[3] = {0, 0, 0};
8745 ParallelIo::size_type bcOffset[3] = {0, 0, 0};
8747 bcCells[0] =
m_grid->getMyBlockNoCells(0);
8749 bcCells[2] =
m_grid->getMyBlockNoCells(2);
8751 MInt noCellsBC = bcCells[0] * bcCells[1] * bcCells[2];
8755 stringstream restartFileName;
8760 stringstream pathStr;
8761 pathStr <<
"/block" <<
m_blockId <<
"/bc2601" << endl;
8762 std::string path = pathStr.str();
8764 for(
MInt var = 0; var <
PV->noVariables; var++) {
8770 "tmpRestartVars[0]");
8773 MInt startGC[3] = {0, 0, 0};
8774 MInt endGC[3] = {0, 0, 0};
8789 for(
MInt i = startGC[2]; i <
m_nCells[2] - endGC[2]; i++) {
8791 for(
MInt k = startGC[0]; k <
m_nCells[0] - endGC[0]; k++) {
8799 for(
MInt var = 0; var <
PV->noVariables; var++) {
8813 for(
MInt var = 0; var <
PV->noVariables; var++) {
8831 for(
MInt var = 0; var <
PV->noVariables; var++) {
8847 if(!isPrimitiveOutput &&
domainId() == 0) {
8848 cout <<
"Restart file has conservative variables, converting STG variables to primitive!" << endl;
8851 stringstream restartFileName;
8856 stringstream blockNumber;
8858 MString blockPathStr =
"/block";
8859 blockPathStr += blockNumber.str();
8861 MInt restartNoEddies = 0;
8865 m_log <<
"STG: NRAN in restart file (" << restartNoEddies <<
") not the same as given in property file ("
8866 <<
m_stgMaxNoEddies <<
"), creating new random distribution of eddies!" << endl;
8869 m_log <<
"STG: Reading in " << restartNoEddies <<
" eddies from restart file" << endl;
8880 FQ->loadedFromRestartFile[
FQ->NU_T] =
true;
8884 ParallelIo::size_type VBStart = 0;
8892 cout <<
"NRAN in property differs from NRAN in restart file"
8893 <<
" NOT READING EDDIES FROM RESTART!" << endl;
8896 MString stgGlobalPathStr =
"stgGlobal";
8898 if(pio.
hasDataset(
"FQeddies", stgGlobalPathStr)) {
8900 cout <<
"FQeddies field is at new position /stgGlobal/FQeddies" << endl;
8904 cout <<
"FQeddies field is NOT at new position! Using old path within solver..." << endl;
8906 stgGlobalPathStr = blockPathStr;
8908 cout <<
"Loading FQeddies from path " << stgGlobalPathStr << endl;
8914 cout <<
"Loading STG Eddies..." << endl;
8923 cout <<
"Loading STG Eddies... SUCCESSFUL!" << endl;
8933 ParallelIo::size_type bcActiveCells[3] = {
m_grid->getMyBlockNoCells(0),
m_grid->getMyBlockNoCells(1), 3};
8936 ParallelIo::size_type bcOffset[3] = {0, 0, 0};
8937 const MInt noCellsBC = bcCells[0] * bcCells[1] * bcCells[2];
8941 cout <<
"Loading STG Datasets..." << endl;
8947 stringstream fieldName;
8948 stringstream stgPath;
8949 stgPath << blockPathStr <<
"/stg";
8950 fieldName <<
"stgFQ" << var;
8951 pioLocal.
readArray(&tmpRestartVars[var * noCellsBC], stgPath.str(), fieldName.str(),
nDim, bcOffset,
8956 for(
MInt k = (bcActiveCells[0] - 1); k >= 0; k--) {
8957 for(
MInt j = (bcActiveCells[1] - 1); j >= 0; j--) {
8958 for(
MInt i = (bcActiveCells[2] - 1); i >= 0; i--) {
8959 const MInt cellId_org = i + (j + k * bcActiveCells[1]) * bcActiveCells[2];
8960 const MInt i_new = i;
8963 const MInt cellId = i_new + (j_new + k_new * bcCells[1]) * bcCells[2];
8966 tmpRestartVars[var * noCellsBC + cellId] = tmpRestartVars[var * noCellsBC + cellId_org];
8967 tmpRestartVars[var * noCellsBC + cellId_org] = F0;
8974 for(
MInt j = 0; j < bcCells[1]; j++) {
8975 for(
MInt i = 0; i < bcCells[2]; i++) {
8976 const MInt gcId0 = i + (j + (0) * bcCells[1]) * bcCells[2];
8977 const MInt acId0 = i + (j + (bcCells[0] - 4) * bcCells[1]) * bcCells[2];
8978 const MInt gcId1 = i + (j + (1) * bcCells[1]) * bcCells[2];
8979 const MInt acId1 = i + (j + (bcCells[0] - 3) * bcCells[1]) * bcCells[2];
8981 const MInt gcId2 = i + (j + (bcCells[0] - 2) * bcCells[1]) * bcCells[2];
8982 const MInt acId2 = i + (j + (2) * bcCells[1]) * bcCells[2];
8983 const MInt gcId3 = i + (j + (bcCells[0] - 1) * bcCells[1]) * bcCells[2];
8984 const MInt acId3 = i + (j + (3) * bcCells[1]) * bcCells[2];
8987 tmpRestartVars[var * noCellsBC + gcId0] = tmpRestartVars[var * noCellsBC + acId0];
8988 tmpRestartVars[var * noCellsBC + gcId1] = tmpRestartVars[var * noCellsBC + acId1];
8989 tmpRestartVars[var * noCellsBC + gcId2] = tmpRestartVars[var * noCellsBC + acId2];
8990 tmpRestartVars[var * noCellsBC + gcId3] = tmpRestartVars[var * noCellsBC + acId3];
8996 for(
MInt k = 0; k < bcCells[0]; k++) {
8997 for(
MInt i = 0; i < bcCells[2]; i++) {
8998 const MInt gc1 = i + (1 + k * bcCells[1]) * bcCells[2];
8999 const MInt gc2 = i + (0 + k * bcCells[1]) * bcCells[2];
9000 const MInt ac1 = i + (2 + k * bcCells[1]) * bcCells[2];
9001 const MInt ac2 = i + (3 + k * bcCells[1]) * bcCells[2];
9004 tmpRestartVars[var * noCellsBC + gc1] =
9005 2 * tmpRestartVars[var * noCellsBC + ac1] - tmpRestartVars[var * noCellsBC + ac2];
9006 tmpRestartVars[var * noCellsBC + gc2] =
9007 2 * tmpRestartVars[var * noCellsBC + gc1] - tmpRestartVars[var * noCellsBC + ac1];
9013 for(
MInt k = 0; k < bcCells[0]; k++) {
9014 for(
MInt i = 0; i < bcCells[2]; i++) {
9015 const MInt gc1 = i + (bcCells[1] - 2 + k * bcCells[1]) * bcCells[2];
9016 const MInt gc2 = i + (bcCells[1] - 1 + k * bcCells[1]) * bcCells[2];
9017 const MInt ac1 = i + (bcCells[1] - 3 + k * bcCells[1]) * bcCells[2];
9018 const MInt ac2 = i + (bcCells[1] - 4 + k * bcCells[1]) * bcCells[2];
9022 tmpRestartVars[var * noCellsBC + gc1] =
9023 2 * tmpRestartVars[var * noCellsBC + ac1] - tmpRestartVars[var * noCellsBC + ac2];
9024 tmpRestartVars[var * noCellsBC + gc2] =
9025 2 * tmpRestartVars[var * noCellsBC + gc1] - tmpRestartVars[var * noCellsBC + ac1];
9035 cout <<
"Loading STG Datasets... SUCCESSFUL!" << endl;
9043 for(
MInt i = 0; i < 3; i++) {
9048 MInt cellIdBCGlobal = globalI + (globalJ + globalK * bcCells[1]) * bcCells[2];
9053 m_cells->
stg_fq[var][cellIdBC] = tmpRestartVars[var * noCellsBC + cellIdBCGlobal];
9056 if(!isPrimitiveOutput) {
9064 const MFloat u = rhoU / rho;
9065 const MFloat v = rhoV / rho;
9066 const MFloat w = rhoW / rho;
9067 const MFloat p = gammaMinusOne * (rhoE - F1B2 * rho * (
POW2(u) +
POW2(v) +
POW2(w)));
9115 const MFloat dudxi = u[IPJK] - u[IMJK];
9116 const MFloat dudet = u[IJPK] - u[IJMK];
9117 const MFloat dudze = u[IJKP] - u[IJKM];
9119 const MFloat dvdxi = v[IPJK] - v[IMJK];
9120 const MFloat dvdet = v[IJPK] - v[IJMK];
9121 const MFloat dvdze = v[IJKP] - v[IJKM];
9123 const MFloat dwdxi = w[IPJK] - w[IMJK];
9124 const MFloat dwdet = w[IJPK] - w[IJMK];
9125 const MFloat dwdze = w[IJKP] - w[IJKM];
9148 vortx[IJK] = F1B2 * (dwdy - dvdz) / jac[IJK];
9149 vorty[IJK] = F1B2 * (dudz - dwdx) / jac[IJK];
9150 vortz[IJK] = F1B2 * (dvdx - dudy) / jac[IJK];
9165 MFloat d[3] = {F0, F0, F0};
9166 MFloat e[3] = {F0, F0, F0};
9188 const MFloat dudx = FcellJac
9192 const MFloat dudy = FcellJac
9196 const MFloat dudz = FcellJac
9201 const MFloat dvdx = FcellJac
9205 const MFloat dvdy = FcellJac
9209 const MFloat dvdz = FcellJac
9214 const MFloat dwdx = FcellJac
9218 const MFloat dwdy = FcellJac
9222 const MFloat dwdz = FcellJac
9228 J(0, 0) =
POW2(dudx) + dvdx * dudy + dwdx * dudz;
9229 J(0, 1) = F1B2 * (dudx * dudy + dudx * dvdx + dvdy * dudy + dvdy * dvdx + dwdy * dudz + dvdz * dwdx);
9230 J(0, 2) = F1B2 * (dudx * dudz + dudx * dwdx + dudy * dvdz + dvdx * dwdy + dwdz * dudz + dwdz * dwdx);
9232 J(1, 1) =
POW2(dvdy) + dvdx * dudy + dvdz * dwdy;
9233 J(1, 2) = F1B2 * (dudz * dvdx + dwdx * dudy + dvdy * dvdz + dvdy * dwdy + dwdz * dvdz + dwdz * dwdy);
9236 J(2, 2) =
POW2(dwdz) + dwdx * dudz + dvdz * dwdy;
9246 lambda2[cellId] = d[1];
9265 for(
MInt var = 0; var <
PV->noVariables; var++) {
9277 MPI_SUM,
m_StructuredComm, AT_,
"m_intpPointsHasPartnerLocal",
"m_intpPointsHasPartnerGlobal");
9283 for(
MInt var = 0; var <
PV->noVariables; var++) {
9291 "m_intpPointsVarsLocal[0][0]",
"m_intpPointsVarsGlobal[0][0]");
9296 for(
MInt var = 0; var <
PV->noVariables; var++) {
9310 stringstream fileName;
9321 ParallelIo::size_type dataOffset[2] = {0, 0};
9325 ParallelIo::size_type noDims = 1;
9326 dataSize[0] = dataSize[1];
9333 MString solutionpath =
"field";
9334 solutionpath += path.str();
9335 const char* dsetname = solutionpath.c_str();
9337 for(
MInt v = 0; v <
PV->noVariables; v++) {
9346 for(
MInt v = 0; v <
PV->noVariables; v++) {
9358 for(
MInt v = 0; v <
PV->noVariables; v++) {
9362 pio.
writeArray(&empty, dsetname,
"x", noDims, dataOffset, dataSize);
9363 pio.
writeArray(&empty, dsetname,
"y", noDims, dataOffset, dataSize);
9364 pio.
writeArray(&empty, dsetname,
"z", noDims, dataOffset, dataSize);
9399 stringstream blockNumber;
9401 MString blockPathStr =
"/block";
9402 blockPathStr += blockNumber.str();
9406 for(
MInt var = 0; var <
PV->noVariables; var++) {
9417 const MInt IPJK = IJK + INC[0];
9418 const MInt IMJK = IJK - INC[0];
9419 const MInt IJPK = IJK + INC[1];
9420 const MInt IJMK = IJK - INC[1];
9421 const MInt IJKP = IJK + INC[2];
9422 const MInt IJKM = IJK - INC[2];
9424 const MFloat dvardxi = F1B2 * (var[IPJK] - var[IMJK]);
9425 const MFloat dvardet = F1B2 * (var[IJPK] - var[IJMK]);
9426 const MFloat dvardze = F1B2 * (var[IJKP] - var[IJKM]);
9438 const MInt IPJK = IJK + INC[0];
9439 const MInt IMJK = IJK - INC[0];
9440 const MInt IJPK = IJK + INC[1];
9441 const MInt IJMK = IJK - INC[1];
9442 const MInt IJKP = IJK + INC[2];
9443 const MInt IJKM = IJK - INC[2];
9445 const MFloat dvardxi = var[IPJK] - var[IMJK];
9446 const MFloat dvardet = var[IJPK] - var[IJMK];
9447 const MFloat dvardze = var[IJKP] - var[IJKM];
9467 m_log <<
"loading average restart file ... " << endl;
9469 cout <<
"Opening file: " << fileName << endl;
9473 MString fileNameStr = fileName;
9477 m_log <<
"-> reading in the data ... " << endl;
9479 m_log <<
"Current number of postprocessing samples: " <<
m_noSamples << endl;
9480 stringstream blockNumber;
9482 MString blockPathStr =
"/block";
9483 blockPathStr += blockNumber.str();
9489 pio.
readArray(sum[0], blockPathStr,
"u",
nDim, ioOffset, ioSize);
9490 pio.
readArray(sum[1], blockPathStr,
"v",
nDim, ioOffset, ioSize);
9491 pio.
readArray(sum[2], blockPathStr,
"w",
nDim, ioOffset, ioSize);
9492 pio.
readArray(sum[3], blockPathStr,
"rho",
nDim, ioOffset, ioSize);
9493 pio.
readArray(sum[4], blockPathStr,
"p",
nDim, ioOffset, ioSize);
9505 pio.
readArray(sum[offset + 0], blockPathStr,
"vortx",
nDim, ioOffset, ioSize);
9506 pio.
readArray(sum[offset + 1], blockPathStr,
"vorty",
nDim, ioOffset, ioSize);
9507 pio.
readArray(sum[offset + 2], blockPathStr,
"vortz",
nDim, ioOffset, ioSize);
9510 pio.
readArray(square[0], blockPathStr,
"uu",
nDim, ioOffset, ioSize);
9511 pio.
readArray(square[1], blockPathStr,
"vv",
nDim, ioOffset, ioSize);
9512 pio.
readArray(square[2], blockPathStr,
"ww",
nDim, ioOffset, ioSize);
9513 pio.
readArray(square[3], blockPathStr,
"uv",
nDim, ioOffset, ioSize);
9514 pio.
readArray(square[4], blockPathStr,
"vw",
nDim, ioOffset, ioSize);
9515 pio.
readArray(square[5], blockPathStr,
"uw",
nDim, ioOffset, ioSize);
9517 pio.
readArray(square[6], blockPathStr,
"pp",
nDim, ioOffset, ioSize);
9520 pio.
readArray(square[7], blockPathStr,
"vortxvortx",
nDim, ioOffset, ioSize);
9521 pio.
readArray(square[8], blockPathStr,
"vortyvorty",
nDim, ioOffset, ioSize);
9522 pio.
readArray(square[9], blockPathStr,
"vortzvortz",
nDim, ioOffset, ioSize);
9526 pio.
readArray(cube[0], blockPathStr,
"uuu",
nDim, ioOffset, ioSize);
9527 pio.
readArray(cube[1], blockPathStr,
"vvv",
nDim, ioOffset, ioSize);
9528 pio.
readArray(cube[2], blockPathStr,
"www",
nDim, ioOffset, ioSize);
9532 pio.
readArray(fourth[0], blockPathStr,
"uuuu",
nDim, ioOffset, ioSize);
9533 pio.
readArray(fourth[1], blockPathStr,
"vvvv",
nDim, ioOffset, ioSize);
9534 pio.
readArray(fourth[2], blockPathStr,
"wwww",
nDim, ioOffset, ioSize);
9537 m_log <<
"loading Restart file ... SUCCESSFUL " << endl;
9548 m_log <<
"loading averaged variables file ... " << endl;
9551 MString fileNameStr = fileName;
9554 stringstream blockNumber;
9556 MString blockPathStr =
"/block";
9557 blockPathStr += blockNumber.str();
9571 m_log <<
"loading Restart file ... SUCCESSFUL " << endl;
9582 MInt cellId_org = 0;
9584 MInt i_new, j_new, k_new;
9604 m_favre[var][cellId_org] = F0;
9614 for(
MInt var = 0; var <
nDim; var++) {
9616 m_cube[var][cellId_org] = F0;
9621 for(
MInt var = 0; var <
nDim; var++) {
9638 MInt cellId_org = 0;
9640 MInt i_new, j_new, k_new;
9660 m_favre[var][cellId_org] = F0;
9668template <MFloat (FvStructuredSolver<3>::*pressure_func)(MInt) const>
9679 const MFloat fRho = F1 / cvars[
CV->RHO][cellId];
9681 for(
MInt vel = 0; vel <
nDim; ++vel) {
9682 pvars[vel][cellId] = cvars[vel][cellId] * fRho;
9683 velPOW2 +=
POW2(pvars[vel][cellId]);
9687 pvars[
PV->RHO][cellId] = cvars[
CV->RHO][cellId];
9688 pvars[
PV->P][cellId] =
9690 * (cvars[
CV->RHO_E][cellId] - F1B2 * pvars[
PV->RHO][cellId] * velPOW2 + (this->*pressure_func)(cellId));
9695 cvars[
CV->RANS_VAR[ransVar]][cellId] =
mMax(cvars[
CV->RANS_VAR[ransVar]][cellId], F0);
9696 pvars[
PV->RANS_VAR[ransVar]][cellId] = cvars[
CV->RANS_VAR[ransVar]][cellId] * fRho;
9705 computePrimitiveVariables_<&FvStructuredSolver::pressure_twoEqRans>();
9725 "ReconstructionConstants", 0.0, AT_);
9736 if(lx % 2 != 0 || ly % 2 != 0 || lz % 2 != 0) {
9737 stringstream errorMessage;
9738 errorMessage <<
" FFTInit: Domain size must NOT be an odd number!: (lx)x(ly)x(lz)-> " << lx <<
"x" << ly <<
"x"
9740 mTerm(1, AT_, errorMessage.str());
9743 m_log <<
" --- initializing FFTW --- " << endl;
9744 m_log <<
" domain size = " << lx <<
"x" << ly <<
"x" << lz << endl;
9746 MFloat waveVector[3], k0;
9748 complex<MFloat>* fourierCoefficient =
new complex<MFloat>[3];
9750 fftw_complex *uHatField, *vHatField, *wHatField;
9752 fftw_plan planU, planV, planW;
9755 uHatField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
9756 vHatField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
9757 wHatField = (fftw_complex*)fftw_malloc(lx * ly * lz *
sizeof(fftw_complex));
9760 planU = fftw_plan_dft_3d(lx, ly, lz, uHatField, uPhysField, FFTW_BACKWARD, FFTW_MEASURE);
9761 planV = fftw_plan_dft_3d(lx, ly, lz, vHatField, vPhysField, FFTW_BACKWARD, FFTW_MEASURE);
9762 planW = fftw_plan_dft_3d(lx, ly, lz, wHatField, wPhysField, FFTW_BACKWARD, FFTW_MEASURE);
9764 for(
MInt p = 0; p < lx; p++) {
9765 for(
MInt q = 0; q < ly; q++) {
9766 for(
MInt r = 0; r < lz; r++) {
9767 MInt cellId = r + lz * (q + ly * p);
9769 uHatField[cellId][0] = 0.0;
9770 uHatField[cellId][1] = 0.0;
9771 uPhysField[cellId][0] = 0.0;
9772 uPhysField[cellId][1] = 0.0;
9775 vHatField[cellId][0] = 0.0;
9776 vHatField[cellId][1] = 0.0;
9777 vPhysField[cellId][0] = 0.0;
9778 vPhysField[cellId][1] = 0.0;
9781 wHatField[cellId][0] = 0.0;
9782 wHatField[cellId][1] = 0.0;
9783 wPhysField[cellId][0] = 0.0;
9784 wPhysField[cellId][1] = 0.0;
9799 k0 = 2.0 * PI / (lx / noPeakModes);
9801 for(
MInt p = 0; p <= lx / 2; p++) {
9802 for(
MInt q = 0; q <= ly / 2; q++) {
9803 for(
MInt r = 0; r <= lz / 2; r++) {
9805 waveVector[0] = (p)*2.0 * PI / lx;
9806 waveVector[1] = (q)*2.0 * PI / ly;
9807 waveVector[2] = (r)*2.0 * PI / lz;
9812 uHatField[r + lz * (q + ly * p)][0] = real(fourierCoefficient[0]);
9813 uHatField[r + lz * (q + ly * p)][1] = imag(fourierCoefficient[0]);
9815 vHatField[r + lz * (q + ly * p)][0] = real(fourierCoefficient[1]);
9816 vHatField[r + lz * (q + ly * p)][1] = imag(fourierCoefficient[1]);
9818 wHatField[r + lz * (q + ly * p)][0] = real(fourierCoefficient[2]);
9819 wHatField[r + lz * (q + ly * p)][1] = imag(fourierCoefficient[2]);
9822 if(p > 1 && q > 1 && r > 1) {
9823 if(p < lx / 2 && q < ly / 2 && r < lz / 2) {
9826 uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = uHatField[r + lz * (q + ly * p)][0];
9827 uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -uHatField[r + lz * (q + ly * p)][1];
9829 vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = vHatField[r + lz * (q + ly * p)][0];
9830 vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -vHatField[r + lz * (q + ly * p)][1];
9832 wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = wHatField[r + lz * (q + ly * p)][0];
9833 wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -wHatField[r + lz * (q + ly * p)][1];
9844 fftw_execute(planU);
9845 fftw_execute(planV);
9846 fftw_execute(planW);
9850 for(
MInt p = 0; p < lx; p++) {
9851 for(
MInt q = 0; q < ly; q++) {
9852 for(
MInt r = 0; r < lz; r++) {
9853 uPhysField[r + lz * (q + ly * p)][0] /= sqrt(
MFloat(lx * ly * lz));
9854 vPhysField[r + lz * (q + ly * p)][0] /= sqrt(
MFloat(lx * ly * lz));
9855 wPhysField[r + lz * (q + ly * p)][0] /= sqrt(
MFloat(lx * ly * lz));
9857 uPhysField[r + lz * (q + ly * p)][1] /= sqrt(
MFloat(lx * ly * lz));
9858 vPhysField[r + lz * (q + ly * p)][1] /= sqrt(
MFloat(lx * ly * lz));
9859 wPhysField[r + lz * (q + ly * p)][1] /= sqrt(
MFloat(lx * ly * lz));
9864 fftw_destroy_plan(planU);
9865 fftw_destroy_plan(planV);
9866 fftw_destroy_plan(planW);
9867 fftw_free(uHatField);
9868 fftw_free(vHatField);
9869 fftw_free(wHatField);
9882 MFloat r[6], s[6], kAbs, energy;
9883 complex<MFloat> uHat, vHat, wHat;
9886 kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
9892 fourierCoefficient[0] = complex<MFloat>(0, 0);
9893 fourierCoefficient[1] = complex<MFloat>(0, 0);
9894 fourierCoefficient[2] = complex<MFloat>(0, 0);
9898 energy = pow(kAbs / k0, 4.0) * exp(-2.0 * (kAbs / k0) * (kAbs / k0));
9899 energy *= exp(2.0) * 0.499 * (
m_Ma * LBCS)
9907 for(
MInt i = 0; i < 6; i++) {
9908 r[i] =
randnormal(0.0, PI * sqrt(energy) / (SQRT2 * kAbs));
9910 s[i] =
randnormal(0.0, PI * sqrt(energy) / (SQRT2 * kAbs));
9914 uHat = (1.0 - k[0] * k[0] / (kAbs * kAbs)) * complex<MFloat>(r[0] + r[3], s[0] - s[3])
9915 - k[0] * k[1] / (kAbs * kAbs) * complex<MFloat>(r[1] + r[4], s[1] - s[4])
9916 - k[0] * k[2] / (kAbs * kAbs) * complex<MFloat>(r[2] + r[5], s[2] - s[5]);
9919 vHat = -k[1] * k[0] / (kAbs * kAbs) * complex<MFloat>(r[0] + r[3], s[0] - s[3])
9920 + (1.0 - k[1] * k[1] / (kAbs * kAbs)) * complex<MFloat>(r[1] + r[4], s[1] - s[4])
9921 - k[1] * k[2] / (kAbs * kAbs) * complex<MFloat>(r[2] + r[5], s[2] - s[5]);
9924 wHat = -k[2] * k[0] / (kAbs * kAbs) * complex<MFloat>(r[0] + r[3], s[0] - s[3])
9925 - k[2] * k[1] / (kAbs * kAbs) * complex<MFloat>(r[1] + r[4], s[1] - s[4])
9926 + (1.0 - k[2] * k[2] / (kAbs * kAbs)) * complex<MFloat>(r[2] + r[5], s[2] - s[5]);
9934 fourierCoefficient[0] = uHat;
9935 fourierCoefficient[1] = vHat;
9936 fourierCoefficient[2] = wHat;
9942 MInt nghbr[30], dim = 0;
9946 MInt maxNoSingularityRecNghbrIds = 14;
9956 MInt totalCells = 1;
9960 for(
MInt n = 0; n <
nDim + 1; ++n) {
9968 if(len1[j] != 0) totalCells *= len1[j];
9982 for(
MInt kk = start[2]; kk < end[2]; ++kk) {
9983 for(
MInt jj = start[1]; jj < end[1]; ++jj) {
9984 for(
MInt ii = start[0]; ii < end[0]; ++ii) {
9990 nghbr[count++] =
cellIndex(ii + temp[0], jj + temp[1], kk + temp[2]);
9999 nghbr[count++] =
cellIndex(ii + change[0], jj + change[1], kk + change[2]);
10000 nghbr[count++] =
cellIndex(ii + temp[0] + change[0], jj + temp[1] + change[1], kk + temp[2] + change[2]);
10004 cerr <<
"Something wrong with the singularities in the LS coeffiecient computation" << endl;
10011 for(
MInt n = 0; n < count; n++) {
10012 MInt nghbrId = nghbr[n];
10018 weights[n] = 1 / dxdx;
10021 MInt id2 = ii - start[0] + ((jj - start[1]) + (kk - start[2]) * len1[1]) * len1[0];
10026 maxc =
mMax(maxc, condNum);
10028 if(condNum < F0 || condNum > 1e7 || std::isnan(condNum)) {
10029 cerr <<
domainId() <<
" SVD decomposition for pointId " << ijk
10030 <<
" with large condition number: " << condNum <<
" num of neighbor" << count <<
"x" << recDim <<
" "
10031 <<
" coords " <<
m_grid->m_coordinates[0][ijk] <<
", " <<
m_grid->m_coordinates[1][ijk] <<
", "
10032 <<
m_grid->m_coordinates[2][ijk] << endl;
10044 mTerm(1,
"This 3D version is not tested yet!");
10065 const MInt pp[3][12] = {{0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1},
10066 {0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1},
10067 {0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0}};
10071 const MInt normalDir = face / 2;
10072 const MInt firstTangentialDir = (normalDir + 1) %
nDim;
10073 const MInt secondTangentialDir = (normalDir + 2) %
nDim;
10074 const MInt normalDirStart = start[normalDir];
10075 const MInt firstTangentialStart = start[firstTangentialDir];
10076 const MInt firstTangentialEnd = end[firstTangentialDir];
10077 const MInt secondTangentialStart = start[secondTangentialDir];
10078 const MInt secondTangentialEnd = end[secondTangentialDir];
10079 const MInt incp[
nDim] = {IJKP[normalDir], IJKP[firstTangentialDir], IJKP[secondTangentialDir]};
10080 const MInt inc[
nDim] = {IJ[normalDir], IJ[firstTangentialDir], IJ[secondTangentialDir]};
10082 const MInt n = (face % 2) * 2 - 1;
10083 const MInt g1p = normalDirStart + 2 * ((
MInt)(0.5 - (0.5 * (
MFloat)n)));
10084 const MInt g1 = normalDirStart + (
MInt)(0.5 - (0.5 * (
MFloat)n));
10086 for(
MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
10087 for(
MInt t2 = secondTangentialStart; t2 < secondTangentialEnd; t2++) {
10088 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1] + t2 * inc[2];
10091 const MInt ijk = g1p * incp[0] + t1 * incp[1] + t2 * incp[2];
10096 helper[normalDir] = 1;
10105 for(
MInt dim = 0; dim <
nDim; dim++) {
10106 firstVec[dim] =
m_grid->m_coordinates[dim][pp2] -
m_grid->m_coordinates[dim][pp1];
10107 secondVec[dim] =
m_grid->m_coordinates[dim][pp3] -
m_grid->m_coordinates[dim][pp1];
10108 normalVec_[dim] =
m_grid->m_coordinates[dim][pp3_] -
m_grid->m_coordinates[dim][pp1];
10113 const MFloat normalLength = sqrt(
POW2(normalVec[0]) +
POW2(normalVec[1]) +
POW2(normalVec[2]));
10115 MFloat sgn = (std::inner_product(&normalVec[0], &normalVec[0] +
nDim, &normalVec_[0], 0.0) < 0.0) ? -1 : 1;
10118 for(
MInt dim = 0; dim <
nDim; dim++) {
10119 normalVec[dim] /= normalLength;
10120 m_cells->
fq[
FQ->NORMAL[dim]][cellIdG1] = sgn * normalVec[dim];
10129 mTerm(1,
"Not tested yet!");
10132 if(singularity.BC == -6000) {
10133 MBool takeIt =
false;
10134 for(
MInt n = 0; n < singularity.Nstar; ++n) {
10135 if(singularity.BCsingular[n] == -6002) takeIt =
true;
10141 if(singularity.end[n] - singularity.start[n] > 1) {
10144 start[n] = singularity.
start[n] + 1;
10145 end[n] = singularity.end[n] - 1;
10147 start[n] = singularity.start[n];
10148 end[n] = singularity.end[n];
10152 for(
MInt kk = start[2]; kk < end[2]; ++kk) {
10153 for(
MInt jj = start[1]; jj < end[1]; ++jj) {
10154 for(
MInt ii = start[0]; ii < end[0]; ++ii) {
10157 for(
MInt m = 0; m < 2; ++m) {
10158 const MInt* change = singularity.displacement[m];
10159 const MInt nghbr =
cellIndex(ii + change[0], jj + change[1], kk + change[2]);
10176 std::vector<MInt> receiveSizes;
10179 if(snd->bcId == -6000 || snd->bcId == -6002) {
10183 for(
MInt dim = 0; dim <
nDim; ++dim)
10184 size *= snd->endInfoCells[dim] - snd->startInfoCells[dim];
10189 for(
MInt k = snd->startInfoCells[2]; k < snd->endInfoCells[2]; k++) {
10190 for(
MInt j = snd->startInfoCells[1]; j < snd->endInfoCells[1]; j++) {
10191 for(
MInt i = snd->startInfoCells[0]; i < snd->endInfoCells[0]; i++) {
10197 snd->cellBuffer[pos] =
m_cells->
fq[
FQ->POROSITY][cellId];
10199 snd->cellBuffer[(1 + d) * size + pos] =
m_cells->
fq[
FQ->NORMAL[d]][cellId];
10208 MInt err =
MPI_Isend((
void*)&snd->cellBuffer[0], size * (
nDim + 1), MPI_DOUBLE, snd->nghbrId, tag,
10210 if(err) cout <<
"rank " <<
domainId() <<
" sending throws error " << endl;
10216 for(
MInt dim = 0; dim <
nDim; ++dim)
10217 size *= snd->endInfoCells[dim] - snd->startInfoCells[dim];
10218 receiveSizes.push_back(size);
10225 std::vector<MFloat> bufferRcv(std::accumulate(receiveSizes.begin(), receiveSizes.end(), 0) * (
nDim + 1));
10230 if(rcv->bcId == -6002 || rcv->bcId == -6000) {
10231 const MInt rcvSize = receiveSizes[cnt];
10234 &rcv->mpi_request, AT_,
"(void*)&rcvSize");
10235 if(err) cout <<
"rank " <<
domainId() <<
" sending throws error " << endl;
10237 offset += rcvSize * (1 +
nDim);
10246 if(snd->bcId == -6002 || snd->bcId == -6000) {
10247 MPI_Wait(&(snd->mpi_request), &(snd->mpi_status), AT_);
10253 if(rcv->bcId == -6002 || rcv->bcId == -6000) {
10254 MPI_Wait(&(rcv->mpi_request), &(rcv->mpi_status), AT_);
10266 if(rcv->bcId == -6002 || rcv->bcId == -6000) {
10269 MInt k2, j2, i2, id2;
10275 MInt totalCells = 1;
10279 len1[j] = rcv->endInfoCells[j] - rcv->startInfoCells[j];
10280 if(len1[j] != 0) totalCells *= len1[j];
10282 step2[rcv->orderInfo[j]] = rcv->stepInfo[j];
10286 ASSERT(totalCells == receiveSizes[cnt],
"");
10290 end2[j] = len1[j] - 1;
10291 len2[rcv->orderInfo[j]] = len1[j];
10294 start2[j] = end2[j];
10301 for(
MInt k = rcv->startInfoCells[2]; k < rcv->endInfoCells[2]; k++) {
10303 for(
MInt j = rcv->startInfoCells[1]; j < rcv->endInfoCells[1]; j++) {
10305 for(
MInt i = rcv->startInfoCells[0]; i < rcv->endInfoCells[0]; i++) {
10306 start1[rcv->orderInfo[0]] = i2;
10307 start1[rcv->orderInfo[1]] = j2;
10308 start1[rcv->orderInfo[2]] = k2;
10310 id2 = start1[0] + (start1[1] + start1[2] * len2[1]) * len2[0];
10312 m_cells->
fq[
FQ->POROSITY][cellId] = bufferRcv[offset * (
nDim + 1) + id2];
10314 normals_temp[
m_noCells * d + cellId] = bufferRcv[offset * (
nDim + 1) + (d + 1) * totalCells + id2];
10326 offset += totalCells;
10335 if(singularity.BC == -6000) {
10336 MBool takeIt =
false;
10337 for(
MInt n = 0; n < singularity.Nstar; ++n) {
10338 if(singularity.BCsingular[n] == -6002) takeIt =
true;
10344 if(singularity.end[n] - singularity.start[n] > 1) {
10347 start[n] = singularity.
start[n] + 1;
10348 end[n] = singularity.end[n] - 1;
10350 start[n] = singularity.start[n];
10351 end[n] = singularity.end[n];
10355 const MInt nstar = singularity.Nstar;
10357 for(
MInt kk = start[2]; kk < end[2]; ++kk) {
10358 for(
MInt jj = start[1]; jj < end[1]; ++jj) {
10359 for(
MInt ii = start[0]; ii < end[0]; ++ii) {
10365 for(
MInt m = 0; m < nstar - 1; ++m) {
10366 const MInt* change = singularity.displacement[m];
10367 const MInt nghbr =
cellIndex(ii + change[0], jj + change[1], kk + change[2]);
10369 temp[d] += normals_temp[
m_noCells * d + nghbr];
10374 m_cells->
fq[
FQ->NORMAL[d]][IJ] = temp[d] / (2 * nstar);
10382 for(
MInt m = 0; m < nstar - 1; ++m) {
10383 const MInt* change = singularity.displacement[m];
10384 const MInt nghbr =
cellIndex(ii + change[0], jj + change[1], kk + change[2]);
10398 std::vector<MFloat> sendBuffer;
10399 std::vector<MInt> sendcounts;
10400 std::vector<MInt> snghbrs;
10401 std::vector<MInt> tags;
10403 if (snd->bcId==-6002) {
10404 snghbrs.push_back(snd->nghbrId);
10405 tags.push_back(
domainId()+(snd->tagHelper)*m_solver->noDomains());
10408 for (
MInt dim = 0; dim <
nDim; ++dim)
10409 size *= snd->endInfoCells[dim] - snd->startInfoCells[dim];
10410 sendcounts.push_back(size);
10412 sendBuffer.resize(pos+size);
10413 for(
MInt j=snd->startInfoCells[1]; j<snd->endInfoCells[1]; j++) {
10414 for(
MInt i=snd->startInfoCells[0]; i<snd->endInfoCells[0]; i++) {
10420 sendBuffer[pos++] =
m_cells->
fq[
FQ->POROSITY][cellId];
10426 const MInt noNeighborDomains = snghbrs.size();
10429 std::vector<MInt> recvcounts(noNeighborDomains);
10430 std::vector<MInt> rdispls(noNeighborDomains);
10438 m_solver->domainId(),
10443 std::vector<MInt> sendcounts2(noNeighborDomains, 1);
10451 m_solver->domainId(),
10456 if (rcv->bcId==-6002) {
10457 const MInt tag = rcv->nghbrId+rcv->tagHelper*m_solver->noDomains();
10459 for (n = 0; n < noNeighborDomains; ++n) {
10460 if (tag==recvTags[n])
10463 if (n==noNeighborDomains)
mTerm(1,
"n == noNeighborDomains");
10465 const MFloat*
const recvBuffer_ = &recvBuffer[rdispls[n]];
10466 const MInt noReceivedElements = recvcounts[n];
10479 len1[j]=rcv->endInfoCells[j] - rcv->startInfoCells[j];
10480 if(len1[j]!=0) totalCells*=len1[j];
10482 step2[rcv->orderInfo[j]]=rcv->stepInfo[j];
10486 ASSERT(noReceivedElements==totalCells,
"noReceivedElements == totalCells");
10491 len2[rcv->orderInfo[j]]=len1[j];
10502 for(
MInt j=rcv->startInfoCells[1]; j<rcv->endInfoCells[1]; j++) {
10504 for(
MInt i=rcv->startInfoCells[0]; i<rcv->endInfoCells[0]; i++) {
10505 start1[rcv->orderInfo[0]]=i2;
10506 start1[rcv->orderInfo[1]]=j2;
10508 id2=start1[0]+start1[1]*len2[0];
10510 m_cells->
fq[
FQ->POROSITY][cellId]= recvBuffer_[id2];
10527 const MInt IP1 = I + IJK[dim];
10528 const MInt IM1 = I - IJK[dim];
10529 const MInt IP2 = I + 2 * IJK[dim];
10538 *
mMax(
mMax(fabs((PIM2 - PIM1) /
mMin(PIM2, PIM1)), fabs((PIM1 - PIP1) /
mMin(PIM1, PIP1))),
10539 fabs((PIP1 - PIP2) /
mMin(PIP1, PIP2))));
MLong allocatedBytes()
Return the number of allocated bytes.
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
static MInt propertyLength(const MString &name, MInt solverId=m_noSolvers)
Returns the number of elements of a property.
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
void computeLambda2Criterion() override
Function to compute the lambda_2 criterion.
void shiftAverageCellValuesRestart()
static constexpr MInt xsd
void viscousFluxCorrection()
MFloat computeTotalPressure()
void initInterpolatedPoints()
Compute all the interpolation coefficients necessary for the interpolation output.
MFloat getCellLengthY(MInt, MInt, MInt) override
void AusmLES_PTHRC(MFloat *QLeft, MFloat *QRight, MInt dim, MInt cellId)
AUSM PTHRC.
MBool maxResidual()
Computation of the maximum residual.
void loadSampleFile(MString) override
Loads primitive variables from an HDF5 file.
void allocateSingularities()
void scatter(const MBool, std::vector< std::unique_ptr< StructuredComm< 3 > > > &) override
void initBodyForce() override
void spanwiseWaveReorder() override
void computePrimitiveVariables() override
void zonalExchange() override
void spanwiseAvgZonal(std::vector< MFloat * > &) override
void manualInterpolationCorrection()
Manually correct errors made by the restart interpolation.
void applyBoundaryCondition() override
void MusclVenkatakrishnan3D()
MUSCL with Venkatakrishan limiter.
void shiftAverageCellValues()
Shifts the averaged variables.
void computeTimeStep() override
MFloat dvardxyz(MInt, MInt, MFloat *) override
void applySandpaperTripAirfoil()
void tripForceCoefficients(MFloat *, MFloat *, MFloat *, MInt, MInt)
MFloat dist(MFloat *a, MFloat *b)
void BARTH_JESPERSON_FCT(MFloat effNghbrDelta, MFloat srfcDelta, MFloat dxEpsSqr, MInt cellPos, MInt var, MFloatScratchSpace &minPhi)
Barth-Jesperson Limiter.
void MusclMinModLimiter()
void gather(const MBool, std::vector< std::unique_ptr< StructuredComm< 3 > > > &) override
static constexpr MInt zsd
void(FvStructuredSolver3D::* reconstructSurfaceData)()
void(FvStructuredSolver3D::* Venkatakrishnan_function)(MFloat, MFloat, MFloat, MInt, MInt, MFloatScratchSpace &)
MInt getCellIdfromCell(const MInt, const MInt, const MInt, const MInt)
void saveNodalBoxes() override
void(FvStructuredSolver3D::* viscFluxMethod)()
MFloat getSampleVorticity(MInt, MInt) override
static constexpr MInt ysd
void interpolateFromDonor()
Interpolates the flow field from a given donor file.
std::vector< std::unique_ptr< StructuredZonalBC > > m_zonalBC
void distributeFluxToCells()
void getFourierCoefficients(MFloat *k, MFloat k0, std::complex< MFloat > *fourierCoefficient)
Generates a single complex coefficient of Fourier series.
void savePointsToAsciiFile(MBool) override
Saves variables of given cells to ASCII file.
void initPointsToAsciiFile() override
void waveSend(std::vector< MPI_Request > &)
void gcExtrapolate(std::vector< MFloat * > &)
void initialCondition() override
Computation of infinity values for the conservative and primitive variables.
MInt getPointIdFromCell(const MInt, const MInt, const MInt)
void getSampleVariables(MInt, MFloat *) override
void VENKATAKRISHNAN_MOD_FCT(MFloat effNghbrDelta, MFloat srfcDelta, MFloat dxEpsSqr, MInt cellPos, MInt var, MFloatScratchSpace &minPhi)
Venkatakrishnan limiter, modified for better results.
MInt cellIndex(const MInt, const MInt, const MInt)
void extrapolateGhostPointCoordinatesBC()
void initSolutionStep(MInt mode)
void applyViscousBoundaryCondition()
void tripFourierCoefficients(MFloat *, MInt, MFloat, MFloat)
MFloat randnormal(MFloat mu, MFloat sigma)
Returns a normal distributed random-number with mu=mean and sigma=standard deviation.
void loadRestartSTG(MBool)
MFloat dvardx(MInt, MFloat *) override
void addGhostPointCoordinateValues()
Extrapolates and exchanges ghost point coordinates.
MFloat computeTotalKineticEngergy()
void computeReconstructionConstantsSVD()
void gcFillGhostCells(std::vector< MFloat * > &)
void Muscl_AusmLES_PTHRC()
void loadAverageRestartFile(const MChar *, MFloat **, MFloat **, MFloat **, MFloat **) override
Loads the postprocessing restart file to continue the postprocessing.
void computePrimitiveVariables_()
void viscousFluxLES()
Viscous flux computation.
FvStructuredSolver3D(MInt, StructuredGrid< 3 > *, MBool *, const MPI_Comm comm)
void Muscl(MInt timerId=-1) override
MInt getPointIdfromPoint(const MInt, const MInt, const MInt, const MInt)
MInt pointIndex(const MInt, const MInt, const MInt)
void computeCumulativeAverage(MBool forceReset) override
void loadAveragedVariables(const MChar *) override
Loads the averaged variables again to do further postprocessing.
void saveInterpolatedPoints() override
void initMovingGrid() override
void AusmDV(MFloat *QLeft, MFloat *QRight, MInt dim, MInt cellId)
void waveReceive(std::vector< MPI_Request > &)
void waveExchange() override
std::unique_ptr< StructuredInterpolation< 3 > > m_structuredInterpolation
void applyInviscidBoundaryCondition()
MFloat getPSI(MInt, MInt)
virtual void computePorousRHS(MBool) override
static constexpr const MInt nDim
void applyBodyForce(const MBool isRestart, const MBool zeroPos) override
std::unique_ptr< FvStructuredSolver3DRans > m_ransSolver
void computeVorticity() override
std::unique_ptr< StructuredBndryCnd< 3 > > m_structuredBndryCnd
void applySandpaperTrip()
void computeZonalConnections()
void moveGrid(const MBool isRestart, const MBool zeroPos) override
void initFFTW(fftw_complex *uPhysField, fftw_complex *vPhysField, fftw_complex *wPhysField, MInt lx, MInt ly, MInt lz, MInt noPeakModes)
void VENKATAKRISHNAN_FCT(MFloat effNghbrDelta, MFloat srfcDelta, MFloat dxEpsSqr, MInt cellPos, MInt var, MFloatScratchSpace &minPhi)
Standard Venkatakrishnan limiter.
void crossProduct(MFloat *, const MFloat *, const MFloat *)
std::unique_ptr< StructuredInterpolation< 3 > > m_pointInterpolation
void computeDomainWidth() override
MFloat getCellCoordinate(MInt, MInt) override
void AusmLES(MFloat *QLeft, MFloat *QRight, const MInt dim, const MInt cellId)
AUSM central.
Base class of the structured solver.
MFloat m_waveEndTransition
MFloat m_wavePressureOutEndTransition
MInt m_pointsToAsciiVarId
MFloat m_waveGradientSuction
MBool m_stgInitialStartup
MFloat ** m_nodalBoxVariables
MFloat m_firstAvrgResidual
virtual MFloat getBlasiusEta(MFloat coordX, MFloat coordY)
MBool isPeriodicComm(std::unique_ptr< StructuredComm< nDim > > &)
std::unique_ptr< FvStructuredSolverWindowInfo< nDim > > m_windowInfo
std::unique_ptr< StructuredFQVariables > FQ
void checkNans()
Checks whole domain for NaNs and adds the number of NaNs globally.
MInt * m_intpPointsOffsets
MInt m_intpPointsNoPointsTotal
MFloat m_waveYBeginTransition
MBool m_auxDataCoordinateLimits
MFloat ** m_pointsToAsciiCoordinates
MFloat m_wavePressureBeginTransition
MInt * m_pointsToAsciiHasPartnerGlobal
MInt m_movingGridStepOffset
MInt m_waveCellsPerWaveLength
MBool m_streamwiseTravelingWave
MFloat m_convergenceCriterion
MBool m_mgExchangeCoordinates
MInt ** m_nodalBoxLocalOffset
void saveSolverSolution(MBool=false, const MBool=false) override
MInt * m_intpPointsNoPoints2D
MInt * m_bc2600noOffsetCells
MInt * m_bc2600noActiveCells
MFloat m_waveBeginTransition
MString m_nodalBoxOutputDir
void loadRestartFile()
Load Restart File (primitive and conservative output) general formulation.
MFloat * m_airfoilNormalVec
MString * m_variableNames
MInt m_pointsToAsciiNoPoints
MFloat * m_waveForceField
MFloat computeRecConstSVD(const MInt ijk, const MInt noNghbrIds, MInt *nghbr, MInt ID, MInt sID, MFloatScratchSpace &tmpA, MFloatScratchSpace &tmpC, MFloatScratchSpace &weights, const MInt recDim)
AUX DATA ENDS /////////////////////////////////////////////////////////////.
MBool m_savePartitionOutput
std::array< MInt, nDim > beginP0()
MBool skipPeriodicDirection(std::unique_ptr< StructuredComm< nDim > > &)
MFloat m_physicalTimeStep
MInt m_nodalBoxTotalLocalSize
MFloat m_waveAmplitudeSuction
void tred2(MFloatScratchSpace &A, MInt dim, MFloat *diag, MFloat *offdiag)
Householder Reduction according to Numercial Recipies in C: The Art of Scientific Computing.
MInt m_waveTimeStepComputed
MFloat m_wavePressureOutBeginTransition
MFloat * m_tripMaxAmpSteady
MInt * m_nodalBoxLocalSize
MFloat m_waveTemporalTransition
std::array< MInt, Timers::_count > m_timers
MInt m_pointsToAsciiOutputInterval
MFloat m_channelInflowPlaneCoordinate
MFloat ** m_bc2600Variables
std::vector< std::unique_ptr< StructuredComm< nDim > > > m_sndComm
StructuredGrid< nDim > * m_grid
MFloat m_wavePressureEndTransition
MInt ** m_nodalBoxLocalDomainOffset
virtual void computeConservativeVariables()
void fixTimeStepTravelingWave()
MString m_waveForceFieldFile
void shiftCellValuesRestart(MBool)
MFloat m_firstMaxResidual
MInt noVariables() const override
std::vector< std::unique_ptr< StructuredComm< nDim > > > m_waveRcvComm
std::array< MInt, nDim > beginP2()
SingularInformation * m_singularity
std::array< MInt, nDim > endM1()
MInt m_outputIterationNumber
std::vector< std::unique_ptr< StructuredComm< nDim > > > m_waveSndComm
MInt m_intpPointsOutputInterval
std::unique_ptr< StructuredInterpolation< nDim > > m_nodalBoxInterpolation
virtual void getBlasiusVelocity(MInt cellId, MFloat *const vel)
Load variables for the specified timeStep.
MFloat ** m_intpPointsVarsLocal
MFloat * m_airfoilWallDist
MFloat m_channelPresOutlet
MString * m_pvariableNames
MFloat m_waveOutEndTransition
MInt ** m_nodalBoxLocalPoints
MInt m_stgNoEddieProperties
MFloat m_wavePenetrationHeight
virtual MFloat getFscEta(MFloat coordX, MFloat coordY)
MBool m_stgCreateNewEddies
std::unique_ptr< MPrimitiveVariables< nDim > > PV
MFloat ** m_pointsToAsciiVars
MFloat m_movingGridTimeOffset
MInt m_pointsToAsciiCounter
MFloat m_zonalAveragingFactor
MInt * m_intpPointsHasPartnerGlobal
std::unique_ptr< MConservativeVariables< nDim > > CV
MBool m_restartInterpolation
MInt * m_nodalBoxTotalLocalOffset
MPI_Comm m_StructuredComm
virtual MFloat getFscPressure(MInt cellId)
MString m_intpPointsOutputDir
void tqli2(MFloat *diag, MFloat *offdiag, MInt dim)
Compute Eigenvalues with implicit shift according to Numercial Recipies in C: The Art of Scientific C...
MInt * m_intpPointsNoPoints
std::array< MInt, nDim > beginP1()
MInt * m_pointsToAsciiHasPartnerLocal
MInt m_pointsToAsciiLastComputationStep
virtual void writePropertiesAsAttributes(ParallelIoHdf5 *pio, MString path)
Overloaded version of writePropertiesAsAttributes that receives ParallelIoHdf5 object pointer instead...
MFloat m_waveYEndTransition
MFloat m_waveGradientPressure
std::unique_ptr< StructuredInterpolation< nDim > > m_pointsToAsciiInterpolation
MBool m_nodalBoxInitialized
void initializeFQField()
Counts the number of necessary FQ fields, allocates them and corrects the indexes of the FQ variable ...
MFloat m_waveAmplitudePressure
void exchange()
SVD STUFF ENDS /////////////////////////////////////////////////////////////.
MInt m_waveNoStepsPerCell
MFloat m_waveOutBeginTransition
MFloat * m_tripMaxAmpFluc
MFloat m_channelPresInlet
MInt * m_intpPointsHasPartnerLocal
MInt m_orderOfReconstruction
std::array< MInt, nDim > endM2()
MFloat ** m_nodalBoxCoordinates
std::vector< std::unique_ptr< StructuredComm< nDim > > > m_rcvComm
MBool m_bc2601InitialStartup
MBool m_nodalBoxWriteCoordinates
MFloat ** m_intpPointsVarsGlobal
MBool m_movingGridInitialStart
MBool m_bc2600InitialStartup
virtual void getFscVelocity(MInt cellId, MFloat *const vel)
Load variables for the specified timeStep.
MFloat ** m_intpPointsCoordinates
MFloat m_waveAmplitudePlus
MInt m_pointsToAsciiLastOutputStep
std::vector< MFloat * > m_zonalSpanwiseAvgVars
MInt m_airfoilNoWallPoints
void allocateAuxDataMaps()
AUX DATA //////////////////////////////////////////////////////////////////.
MInt * m_nodalBoxPartnerLocal
void insertSort(MInt dim, MFloat *list)
Sorting function to sort list in ascending order.
MFloat m_globalDomainWidth
virtual void writeHeaderAttributes(ParallelIoHdf5 *pio, MString fileType)
Overloaded version of writeHeaderAttributes that receives ParallelIoHdf5 object pointer instead of 'f...
void defineArray(maiabd_type type, const MString &name, size_type totalCount)
Create a new array in the file.
void writeArray(const T *array, const MString &name, size_type memoryStride=-1, size_type diskStride=-1)
Write array data to file. [MPI]
MBool hasDataset(const MString &name, MInt dimension)
Check if the file contains an dataset with the given name and dimension.
void getAttribute(T *value, const MString &name, const MString &datasetName="")
Retrieve a file or dataset attribute.
void setAttribute(const T &value, const MString &name, const MString &datasetName="")
Set a file or dataset attribute. [MPI]
void readArray(T *array, const MString &name, size_type memoryStride=-1, size_type diskStride=-1)
Read array data from file. [MPI]
This class is a ScratchSpace.
void fill(T val)
fill the scratch with a given value
MString outputDir() const
Return the directory for output files.
virtual MInt domainId() const
Return the domainId (rank)
MInt m_residualInterval
The number of timesteps before writing the next residual.
const MInt m_solverId
a unique solver identifier
virtual MInt noDomains() const
MFloat m_Re
the Reynolds number
MFloat m_Ma
the Mach number
MFloat ** m_tempWaveSample
MInt getNoPPSquareVars()
Returns number of pp Square variables.
MString * m_avgVariableNames
MString * m_avgFavreNames
MInt getNoPPVars()
Returns number of postprocessing variables.
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
void mTerm(const MInt errorCode, const MString &location, const MString &message)
constexpr Real POW2(const Real x)
MBool approx(const T &, const U &, const T)
constexpr T mMin(const T &x, const T &y)
constexpr T mMax(const T &x, const T &y)
MInt noRansEquations(RansMethod ransMethod)
void printAllocatedMemory(const MLong oldAllocatedBytes, const MString &solverName, const MPI_Comm comm)
Prints currently allocated memory.
MInt globalDomainId()
Return global domain id.
constexpr MLong IPOW2(MInt x)
std::basic_string< char > MString
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Isend
int MPI_Barrier(MPI_Comm comm, const MString &name)
same as MPI_Barrier
int MPI_Allgatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, const int recvcounts[], const int displs[], MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgatherv
int MPI_Irecv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Irecv
int MPI_Wait(MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Wait
int MPI_Waitall(int count, MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Waitall
int MPI_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_Alltoall(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_Alltoall
int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm, const MString &name, const MString &varname)
same as MPI_Bcast
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgather
MFloat distance(const MFloat *a, const MFloat *b)
std::vector< T > mpiExchangePointToPoint(const T *const sendBuffer, const MInt *const snghbrs, const MInt nosnghbrs, const MInt *const sendcounts, const MInt *const rnghbrs, const MInt nornghbrs, const MPI_Comm &mpi_comm, const MInt domainId, const MInt dataBlockSize, MInt *const recvcounts_=nullptr, MInt *const rdispls_=nullptr)
@ MGCellCenterCoordinates