GMGPolar/src/Interpolation/extrapolated_prolongation.cpp at 8e8dc2703974f55c4e470f16e48f69f4b3739761 · SciCompMod/GMGPolar · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#include "../../include/Interpolation/interpolation.h"
using namespace gmgpolar;

/*
 * Extrapolated Prolongation Operator
 * ----------------------------------
 *
 * Extrapolated prolongation is used between the finest most grids in the multigrid hierarchy.
 * It assumes the fine grid comes from a uniform refinement of the coarse grid, so all spacings are equal.
 * Thus fine values can be computed simply by averaging the neighboring coarse nodes.
 *
 * A fine node is classified by the parity of its (i_r, i_theta) indices:
 *
 * 1) (even, even)
 *    Fine node coincides with a coarse node
 *    -> copy value
 *
 * 2) (odd, even)
 *    Node lies between two coarse nodes in radial direction
 *
 *        X ---- O ---- X
 *
 *    -> arithmetic mean of left + right coarse node
 *
 * 3) (even, odd)
 *    Node lies between two coarse nodes in angular direction
 *
 *        X
 *        |
 *        O
 *        |
 *        X
 *
 *    -> arithmetic mean of bottom + top coarse node
 *
 * 4) (odd, odd)
 *    Node lies inside a coarse cell
 *    We extrapolate/average across the diagonal:
 *
 *          X
 *            \
 *              O
 *                \
 *                 X
 *
 *   -> arithmetic mean of bottom right + top left coarse node
 *
 */

static KOKKOS_INLINE_FUNCTION void fineNodeExtrapolatedProlongation(const int i_r, const int i_theta,
                                                                    const PolarGrid<DefaultMemorySpace>& coarse_grid,
                                                                    const PolarGrid<DefaultMemorySpace>& fine_grid,
                                                                    Vector<double>& fine_result,
                                                                    ConstVector<double>& coarse_values)
{
    const int i_r_coarse     = i_r / 2;
    const int i_theta_coarse = i_theta / 2;

    if (i_r & 1) {
        if (i_theta & 1) { /* (odd, odd) -> node in center of coarse cell */
            const double value =
                0.5 * (coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] + /* Bottom right */
                       coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] /* Top left */
                      );
            fine_result[fine_grid.index(i_r, i_theta)] = value;
        }
        else { /* (odd, even) -> between coarse nodes in radial direction */
            const double value = 0.5 * (coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Left */
                                        coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] /* Right */
                                       );
            fine_result[fine_grid.index(i_r, i_theta)] = value;
        }
    }
    else {
        if (i_theta & 1) { /* (even, odd) -> between coarse nodes in angular direction */
            const double value = 0.5 * (coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Bottom */
                                        coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] /* Top */
                                       );
            fine_result[fine_grid.index(i_r, i_theta)] = value;
        }
        else { /* (even, even) -> node lies exactly on coarse grid */
            fine_result[fine_grid.index(i_r, i_theta)] =
                coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)]; /* Center */
        }
    }
}

void Interpolation::applyExtrapolatedProlongation(const PolarGrid<DefaultMemorySpace>& coarse_grid,
                                                  const PolarGrid<DefaultMemorySpace>& fine_grid,
                                                  Vector<double> fine_result, ConstVector<double> coarse_values) const
{
    assert(std::ssize(coarse_values) == coarse_grid.numberOfNodes());
    assert(std::ssize(fine_result) == fine_grid.numberOfNodes());

    /* We split the loops into two regions to better respect the */
    /* access patterns of the smoother and improve cache locality. */

    // The For loop matches circular access pattern */
    Kokkos::parallel_for(
        "Interpolation: Extrapolated Prolongation (Circular)",
        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
            {0, 0}, // Starting point of the index space
            {fine_grid.numberSmootherCircles(), fine_grid.ntheta()} // Ending point of the index space
            ),
        // Kokkos lambda function to execute for each point in the index space
        KOKKOS_LAMBDA(const int i_r, const int i_theta) {
            fineNodeExtrapolatedProlongation(i_r, i_theta, coarse_grid, fine_grid, fine_result, coarse_values);
        });

    /* For loop matches radial access pattern */
    Kokkos::parallel_for(
        "Interpolation: Extrapolated Prolongation (Radial)",
        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
            {0, fine_grid.numberSmootherCircles()}, // Starting point of the index space
            {fine_grid.ntheta(), fine_grid.nr()} // Ending point of the index space
            ),
        // Kokkos lambda function to execute for each point in the index space
        KOKKOS_LAMBDA(const int i_theta, const int i_r) {
            fineNodeExtrapolatedProlongation(i_r, i_theta, coarse_grid, fine_grid, fine_result, coarse_values);
        });

    Kokkos::fence();
}