From 663d1dbbdc47e6d2863eb7b397146b3f41795073 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Wed, 2 Jan 2019 12:41:52 +0100 Subject: Basic grid refinement algorithm by Lagrava et al. Starting point for integration into a more flexible setting. See "Advances in multi-domain lattice Boltzmann grid refinement" (2012) --- apps/adrian/poiseuille2d/poiseuille2d.cpp | 231 ++++++++++++++++++++++++------ 1 file changed, 189 insertions(+), 42 deletions(-) (limited to 'apps/adrian/poiseuille2d') diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp index 1808fbc..8bb982d 100644 --- a/apps/adrian/poiseuille2d/poiseuille2d.cpp +++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp @@ -44,7 +44,7 @@ typedef double T; const T lx = 2.; // length of the channel const T ly = 1.; // height of the channel const int N = 16; // resolution of the model -const T Re = 10.; // Reynolds number +const T Re = 10.; // Reynolds number const T maxPhysT = 40.; // max. simulation time in s, SI unit const T physInterval = 0.25; // interval for the convergence check in s const T residuum = 1e-5; // residuum for the convergence check @@ -274,6 +274,27 @@ public: } }; +template class DESCRIPTOR> +void computeFeq(const Cell& cell, T fEq[DESCRIPTOR::q]) +{ + T rho{}; + T u[2] {}; + cell.computeRhoU(rho, u); + const T uSqr = u[0]*u[0] + u[1]*u[1]; + for (int i=0; i < DESCRIPTOR::q; ++i) { + fEq[i] = lbHelpers::equilibrium(i, rho, u, uSqr); + } +} + +template class DESCRIPTOR> +void computeFneq(const Cell& cell, T fNeq[DESCRIPTOR::q]) +{ + T rho{}; + T u[2] {}; + cell.computeRhoU(rho, u); + lbHelpers::computeFneq(cell, fNeq, rho, u); +} + T order4interpolation(T fm3, T fm1, T f1, T f3) { return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3); @@ -284,8 +305,73 @@ T order3interpolation(T fm1, T f1, T f3) return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3; } +T order2interpolation(T f0, T f1) +{ + return 0.5 * (f0 + f1); +} + +template class DESCRIPTOR> +void storeC2F( + Grid& coarse, + T c2f_rho[], + T c2f_u[][2], + T c2f_fneq[][DESCRIPTOR::q] +) +{ + auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0); + + const int x = coarseLattice.getNx()-2; + for (int y=0; y < coarseLattice.getNy(); ++y) { + T rho{}; + T u[2] {}; + coarseLattice.get(x,y).computeRhoU(rho, u); + c2f_rho[y] = rho; + c2f_u[y][0] = u[0]; + c2f_u[y][1] = u[1]; + + T fNeq[DESCRIPTOR::q] {}; + computeFneq(coarseLattice.get(x,y), fNeq); + for (int i=0; i < DESCRIPTOR::q; ++i) { + c2f_fneq[y][i] = fNeq[i]; + } + } +} + +template class DESCRIPTOR> +void interpolateStoredC2F( + Grid& coarse, + T c2f_rho[], + T c2f_u[][2], + T c2f_fneq[][DESCRIPTOR::q] +) +{ + auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0); + + const int x = coarseLattice.getNx()-2; + for (int y=0; y < coarseLattice.getNy(); ++y) { + T rho{}; + T u[2] {}; + coarseLattice.get(x,y).computeRhoU(rho, u); + c2f_rho[y] = order2interpolation(rho, c2f_rho[y]); + c2f_u[y][0] = order2interpolation(u[0], c2f_u[y][0]); + c2f_u[y][1] = order2interpolation(u[1], c2f_u[y][1]); + + T fNeq[DESCRIPTOR::q] {}; + computeFneq(coarseLattice.get(x,y), fNeq); + for (int i=0; i < DESCRIPTOR::q; ++i) { + c2f_fneq[y][i] = order2interpolation(fNeq[i], c2f_fneq[y][i]); + } + } +} + template class DESCRIPTOR> -void coupleC2F(Grid& coarse, Grid& fine) +void coupleC2F( + Grid& coarse, + Grid& fine, + const T c2f_rho[], + const T c2f_u[][2], + const T c2f_fneq[][DESCRIPTOR::q] +) { auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0); auto& fineLattice = fine.getSuperLattice().getBlockLattice(0); @@ -293,54 +379,104 @@ void coupleC2F(Grid& coarse, Grid& fine) const int x = coarseLattice.getNx()-2; for (int y=0; y < coarseLattice.getNy(); ++y) { + T fEq[DESCRIPTOR::q] {}; + computeFeq(coarseLattice.get(x,y), fEq); + for (int i=0; i < DESCRIPTOR::q; ++i) { - fineLattice.get(0, 2*y)[i] = coarseLattice.get(x,y)[i]; + fineLattice.get(0,2*y)[i] = fEq[i] + coarse.getScalingFactor() * c2f_fneq[y][i]; } } for (int y=2; y < coarseLattice.getNy()-2; ++y) { + const T rho = order4interpolation( + c2f_rho[y-1], + c2f_rho[y], + c2f_rho[y+1], + c2f_rho[y+2] + ); + T u[2] {}; + u[0] = order4interpolation( + c2f_u[y-1][0], + c2f_u[y][0], + c2f_u[y+1][0], + c2f_u[y+2][0] + ); + u[1] = order4interpolation( + c2f_u[y-1][1], + c2f_u[y][1], + c2f_u[y+1][1], + c2f_u[y+2][1] + ); + const T uSqr = u[0]*u[0] + u[1]*u[1]; + for (int i=0; i < DESCRIPTOR::q; ++i) { - fineLattice.get(0,1+2*y)[i] = order4interpolation( - coarseLattice.get(x,y-1)[i], - coarseLattice.get(x,y)[i], - coarseLattice.get(x,y+1)[i], - coarseLattice.get(x,y+2)[i] - ); + const T fneqi = order4interpolation( + c2f_fneq[y-1][i], + c2f_fneq[y][i], + c2f_fneq[y+1][i], + c2f_fneq[y+2][i] + ); + + fineLattice.get(0,1+2*y)[i] = lbHelpers::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi; } } - for (int i=0; i < DESCRIPTOR::q; ++i) { - fineLattice.get(0,1)[i] = order3interpolation( - coarseLattice.get(x,0)[i], - coarseLattice.get(x,1)[i], - coarseLattice.get(x,2)[i] - ); - fineLattice.get(0,fineLattice.getNy()-2)[i] = order3interpolation( - coarseLattice.get(x,coarseLattice.getNy()-1)[i], - coarseLattice.get(x,coarseLattice.getNy()-2)[i], - coarseLattice.get(x,coarseLattice.getNy()-3)[i] - ); + { + const T rho = order3interpolation( + c2f_rho[0], + c2f_rho[1], + c2f_rho[2] + ); + T u[2] {}; + u[0] = order3interpolation( + c2f_u[0][0], + c2f_u[1][0], + c2f_u[2][0] + ); + u[1] = order3interpolation( + c2f_u[0][1], + c2f_u[1][1], + c2f_u[2][1] + ); + const T uSqr = u[0]*u[0] + u[1]*u[1]; + + for (int i=0; i < DESCRIPTOR::q; ++i) { + const T fneqi = order3interpolation( + c2f_fneq[0][i], + c2f_fneq[1][i], + c2f_fneq[2][i] + ); + fineLattice.get(0,1)[i] = lbHelpers::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi; + } } -} -template class DESCRIPTOR> -void computeFneq(const Cell& cell, T fNeq[DESCRIPTOR::q]) -{ - T rho{}; - T u[2]{}; - cell.computeRhoU(rho, u); - lbHelpers::computeFneq(cell, fNeq, rho, u); -} + { + const T rho = order3interpolation( + c2f_rho[coarseLattice.getNy()-1], + c2f_rho[coarseLattice.getNy()-2], + c2f_rho[coarseLattice.getNy()-3] + ); + T u[2] {}; + u[0] = order3interpolation( + c2f_u[coarseLattice.getNy()-1][0], + c2f_u[coarseLattice.getNy()-2][0], + c2f_u[coarseLattice.getNy()-3][0] + ); + u[1] = order3interpolation( + c2f_u[coarseLattice.getNy()-1][1], + c2f_u[coarseLattice.getNy()-2][1], + c2f_u[coarseLattice.getNy()-3][1] + ); + const T uSqr = u[0]*u[0] + u[1]*u[1]; -template class DESCRIPTOR> -void computeFeq(const Cell& cell, T fEq[DESCRIPTOR::q]) -{ - T rho{}; - T u[2]{}; - cell.computeRhoU(rho, u); - const T uSqr = u[0]*u[0] + u[1]*u[1]; - for (int i=0; i < DESCRIPTOR::q; ++i) { - fEq[i] = lbHelpers::equilibrium(i, rho, u, uSqr); + for (int i=0; i < DESCRIPTOR::q; ++i) { + const T fneqi = order3interpolation( + c2f_fneq[coarseLattice.getNy()-1][i], + c2f_fneq[coarseLattice.getNy()-2][i], + c2f_fneq[coarseLattice.getNy()-3][i] + ); + fineLattice.get(0,fineLattice.getNy()-2)[i] = lbHelpers::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi; + } } } @@ -349,7 +485,7 @@ T computeRestrictedFneq(const BlockLatticeView2D& lattice, int x, { T sum = T{0}; for (int i=0; i < DESCRIPTOR::q; ++i) { - T fNeq[DESCRIPTOR::q]{}; + T fNeq[DESCRIPTOR::q] {}; computeFneq(lattice.get(x + DESCRIPTOR::c[i][0], y + DESCRIPTOR::c[i][1]), fNeq); sum += fNeq[i]; } @@ -364,7 +500,7 @@ void coupleF2C(Grid& coarse, Grid& fine) const int x = coarseLattice.getNx()-1; for (int y=0; y < coarseLattice.getNy(); ++y) { - T fEq[DESCRIPTOR::q]{}; + T fEq[DESCRIPTOR::q] {}; computeFeq(fineLattice.get(2,2*y), fEq); const T restrictedFneq = computeRestrictedFneq(fineLattice, 2, 2*y); @@ -443,17 +579,28 @@ int main(int argc, char* argv[]) residuum); timer.start(); + const auto sizeC2F = coarseGrid.getSuperLattice().getBlockLattice(0).getNy(); + T c2f_rho[sizeC2F] {}; + T c2f_u[sizeC2F][2] {}; + T c2f_fneq[sizeC2F][DESCRIPTOR::q] {}; + for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) { if (converge.hasConverged()) { clout << "Simulation converged." << endl; break; } + storeC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq); coarseGrid.getSuperLattice().collideAndStream(); + interpolateStoredC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq); + fineGrid.getSuperLattice().collideAndStream(); - coupleC2F(coarseGrid, fineGrid); + coupleC2F(coarseGrid, fineGrid, c2f_rho, c2f_u, c2f_fneq); + fineGrid.getSuperLattice().collideAndStream(); - coupleC2F(coarseGrid, fineGrid); + storeC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq); + coupleC2F(coarseGrid, fineGrid, c2f_rho, c2f_u, c2f_fneq); + coupleF2C(coarseGrid, fineGrid); getResults( -- cgit v1.2.3