1401 {
1402 TRACE();
1405 cout << "Moving grid" << endl;
1406 }
1408 }
1409
1413
1419
1424
1433
1436
1439 mue[I] = SUTHERLANDLAW(T);
1440 }
1441
1442 uuTilde.fill(F0);
1443 uvTilde.fill(F0);
1444 uwTilde.fill(F0);
1445
1447 cout << "Computing spanwise average/periodic fluctuations" << endl;
1448 }
1449
1452
1453
1457 }
1458 }
1459
1460 vector<MFloat*> spanVariables;
1462 spanVariables.push_back(&spanAvg(var, 0));
1463 }
1464 ppsolver()->spanwiseAvgZonal(spanVariables);
1465 ppsolver()->gcFillGhostCells(spanVariables);
1466
1468 uTilde[I] = u[I] - spanAvg(0, I);
1469 vTilde[I] = v[I] - spanAvg(1, I);
1470 wTilde[I] = w[I] - spanAvg(2, I);
1471 rhoTilde[I] = rho[I] - spanAvg(3, I);
1472 pTilde[I] =
p[I] - spanAvg(4, I);
1473
1474 uuTilde[I] = uTilde[I] * uTilde[I];
1475 uvTilde[I] = uTilde[I] * vTilde[I];
1476 uwTilde[I] = uTilde[I] * wTilde[I];
1477
1478 u[I] = spanAvg(0, I);
1479 v[I] = spanAvg(1, I);
1480 w[I] = spanAvg(2, I);
1481 rho[I] = spanAvg(3, I);
1482 p[I] = spanAvg(4, I);
1483 }
1484 } else {
1485
1486 vector<MFloat*> ppVariables;
1489 }
1490 ppsolver()->spanwiseAvgZonal(ppVariables);
1491 ppsolver()->gcFillGhostCells(ppVariables);
1492 }
1493
1506
1507 dudx.fill(F0);
1508 dudy.fill(F0);
1509 dudz.fill(F0);
1510 dpdx.fill(F0);
1511 duTildedx.fill(F0);
1512 duTildedy.fill(F0);
1513 duTildedz.fill(F0);
1514 dpTildedx.fill(F0);
1515 duuTildedx.fill(F0);
1516 duwTildedz.fill(F0);
1517 duudx.fill(F0);
1518 duwdz.fill(F0);
1519
1524
1525 nududx.fill(F0);
1526 nududz.fill(F0);
1527 nududxx.fill(F0);
1528 nududzz.fill(F0);
1529
1534
1535 nuduTildedx.fill(F0);
1536 nuduTildedz.fill(F0);
1537 nuduTildedxx.fill(F0);
1538 nuduTildedzz.fill(F0);
1539
1541 cout << "Computing gradients" << endl;
1542 }
1543
1545
1547 const MInt* nActiveCells =
ppsolver()->getActiveCells();
1548 const MInt* nOffsetCells =
ppsolver()->getOffsetCells();
1549
1550 for(
MInt i = 1; i < nCells[2] - 1; i++) {
1551 for(
MInt k = 1; k < nCells[0] - 1; k++) {
1552 for(
MInt j = 1; j < nCells[1] - 1; j++) {
1553 const MInt I = i + (j + k * nCells[1]) * nCells[2];
1554 dudx[I] =
ppsolver()->dvardxyz(I, 0, u);
1555 dudy[I] =
ppsolver()->dvardxyz(I, 1, u);
1556 dudz[I] =
ppsolver()->dvardxyz(I, 2, u);
1557 dpdx[I] =
ppsolver()->dvardxyz(I, 0, p);
1558
1559 duTildedx[I] =
ppsolver()->dvardxyz(I, 0, &uTilde[0]);
1560 duTildedy[I] =
ppsolver()->dvardxyz(I, 1, &uTilde[0]);
1561 duTildedz[I] =
ppsolver()->dvardxyz(I, 2, &uTilde[0]);
1562 dpTildedx[I] =
ppsolver()->dvardxyz(I, 0, &pTilde[0]);
1563
1564 duudx[I] =
ppsolver()->dvardxyz(I, 0, uu);
1565 duwdz[I] =
ppsolver()->dvardxyz(I, 2, uw);
1566
1567 duuTildedx[I] =
ppsolver()->dvardxyz(I, 0, &uuTilde[0]);
1568 duwTildedz[I] =
ppsolver()->dvardxyz(I, 2, &uwTilde[0]);
1569
1570 nududx[I] = mue[I] / rho[I] * dudx[I];
1571 nududz[I] = mue[I] / rho[I] * dudz[I];
1572
1573 nuduTildedx[I] = mue[I] / rho[I] * duTildedx[I];
1574 nuduTildedz[I] = mue[I] / rho[I] * duTildedz[I];
1575 }
1576 }
1577 }
1578
1580 cout << "Computing second order gradients" << endl;
1581 }
1582
1583 for(
MInt i = 1; i < nCells[2] - 1; i++) {
1584 for(
MInt k = 1; k < nCells[0] - 1; k++) {
1585 for(
MInt j = 1; j < nCells[1] - 1; j++) {
1586 const MInt I = i + (j + k * nCells[1]) * nCells[2];
1587 nududxx[I] =
ppsolver()->dvardxyz(I, 0, &nududx[0]);
1588 nududzz[I] =
ppsolver()->dvardxyz(I, 2, &nududz[0]);
1589
1590 nuduTildedxx[I] =
ppsolver()->dvardxyz(I, 0, &nuduTildedx[0]);
1591 nuduTildedzz[I] =
ppsolver()->dvardxyz(I, 2, &nuduTildedz[0]);
1592 }
1593 }
1594 }
1595
1596 const MInt globalNoCellsI =
ppsolver()->getGrid()->getMyBlockNoCells(2);
1597 const MInt globalNoCellsK =
ppsolver()->getGrid()->getMyBlockNoCells(0);
1598 const MInt totalNoCellsIK = globalNoCellsI * globalNoCellsK;
1599
1600 const MInt noDecompVars = 11;
1601 MFloatScratchSpace cfDecomposedLocal(noDecompVars, totalNoCellsIK, AT_,
"cfDecomposedLocal");
1602 MFloatScratchSpace cfDecomposedGlobal(noDecompVars, totalNoCellsIK, AT_,
"cfDecomposedGlobal");
1603
1604 cfDecomposedLocal.fill(F0);
1605
1609
1612
1614 cout << "Computing decomposition activeCells[0]: " << nActiveCells[0] << " activeCells[1]: " << nActiveCells[1]
1615 << " activeCells[2]: " << nActiveCells[2] << endl;
1616 cout << "GlobalCellsI: " << globalNoCellsI << endl;
1617 }
1618
1619 for(
MInt i = 0; i < nActiveCells[2]; i++) {
1620 for(
MInt k = 0; k < nActiveCells[0]; k++) {
1621 for(
MInt j = 0; j < nActiveCells[1]; j++) {
1622 const MInt globalId = (i + nOffsetCells[2]) + (k + nOffsetCells[1]) * (globalNoCellsI);
1623 const MInt I = i + noGC + ((j + noGC) + (k + noGC) * nCells[1]) * nCells[2];
1624
1625 const MFloat dy =
ppsolver()->getCellLengthY(i + noGC, j + noGC, k + noGC);
1626
1627 cfDecomposedLocal(0, globalId) += fac * fre0 * mue[I] / rho[I] * dudy[I] * dudy[I] * dy;
1628
1629 cfDecomposedLocal(1, globalId) += fac * fre0 * mue[I] / rho[I] * dudy[I] * duTildedy[I] * dy;
1630
1631 cfDecomposedLocal(2, globalId) += fac * (-uv[I] * dudy[I] * dy);
1632
1633 cfDecomposedLocal(3, globalId) += fac * (-uvTilde[I] * dudy[I] * dy);
1634
1635 cfDecomposedLocal(4, globalId) += fac * ((u[I] - u8) * (u[I] * dudx[I] + v[I] * dudy[I] + w[I] * dudz[I]) * dy);
1636
1637 cfDecomposedLocal(5, globalId) +=
1638 fac * ((u[I] - u8) * (u[I] * duTildedx[I] + v[I] * duTildedy[I] + w[I] * duTildedz[I]) * dy);
1639
1640 cfDecomposedLocal(6, globalId) +=
1641 fac * ((u[I] - u8) * (uTilde[I] * dudx[I] + vTilde[I] * dudy[I] + wTilde[I] * dudz[I]) * dy);
1642
1643 cfDecomposedLocal(7, globalId) += fac * (u[I] - u8) * dpdx[I] * dy;
1644
1645 cfDecomposedLocal(8, globalId) += fac * (u[I] - u8) * dpTildedx[I] * dy;
1646
1647 cfDecomposedLocal(9, globalId) -=
1648 fac * (u[I] - u8) * (fre0 * nududxx[I] + fre0 * nuduTildedxx[I] - duuTildedx[I] - duudx[I]) * dy;
1649 cfDecomposedLocal(10, globalId) -=
1650 fac * (u[I] - u8) * (fre0 * nududzz[I] + fre0 * nuduTildedzz[I] - duwTildedz[I] - duwdz[I]) * dy;
1651 }
1652 }
1653 }
1654
1656 cout << "Before MPI Allreduce" << endl;
1657 }
1658
1660 &cfDecomposedGlobal(0, 0),
1661 noDecompVars * totalNoCellsIK,
1662 MPI_DOUBLE,
1663 MPI_SUM,
1665 AT_,
1666 "cfDecomposedLocal",
1667 "cfDecomposedGlobal");
1668
1669 MFloatScratchSpace cfDecomposedLine(noDecompVars, globalNoCellsI, AT_,
"cfDecomposedLine");
1670 cfDecomposedLine.fill(F0);
1671
1673 cout << "Now averaging in spanwise direction, no cells spanwise: " << globalNoCellsK << endl;
1674 }
1675
1676 for(
MInt var = 0; var < noDecompVars; var++) {
1677 for(
MInt i = 0; i < globalNoCellsI; i++) {
1678 for(
MInt k = 0; k < globalNoCellsK; k++) {
1679 const MInt globalId = i + k * globalNoCellsI;
1680 cfDecomposedLine(var, i) += cfDecomposedGlobal(var, globalId);
1681 }
1682 cfDecomposedLine(var, i) /= (
MFloat)globalNoCellsK;
1683 }
1684 }
1685
1686
1688 MString filename =
"./cf_decomposed.dat";
1689 FILE* f_forces;
1690 f_forces = fopen(filename.c_str(), "w");
1691 for(
MInt i = 0; i < globalNoCellsI; i++) {
1692 for(
MInt var = 0; var < noDecompVars; var++) {
1693 fprintf(f_forces, "%f ", cfDecomposedLine(var, i));
1694 }
1695 fprintf(f_forces, "\n");
1696 }
1697 fclose(f_forces);
1698 }
1699
1701 cout << "Cf decomposition finished!" << endl;
1702 }
1703}
MInt getNoPPVars()
Returns number of postprocessing variables.
MFloat m_sutherlandConstant
MFloat m_sutherlandPlusOne
constexpr Real POW3(const Real x)
constexpr Real POW2(const Real x)
std::basic_string< char > MString
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
constexpr std::underlying_type< FcCell >::type p(const FcCell property)
Converts property name to underlying integer value.
1.9.5