GMGPolar/src/Residual/ResidualGive/residualGive.cpp at 80c320fd68982ebd5672bafe46bb9b309e3734d8 · 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
#include "../../../include/Residual/ResidualGive/residualGive.h"
#include "../../../include/Residual/residual.h"

ResidualGive::ResidualGive(const PolarGrid& grid, const LevelCache& level_cache, const DomainGeometry& domain_geometry,
                           const DensityProfileCoefficients& density_profile_coefficients, bool DirBC_Interior,
                           int num_omp_threads)
    : Residual(grid, level_cache, domain_geometry, density_profile_coefficients, DirBC_Interior, num_omp_threads)
{
}

/* ------------------ */
/* result = rhs - A*x */

// clang-format off
void ResidualGive::computeResidual(Vector<double> result, ConstVector<double> rhs, ConstVector<double> x) const
{
    assert(result.size() == x.size());

    copy_vector(result, rhs);

    /* Single-threaded execution */
    if (num_omp_threads_ == 1) {
        for (int i_r = 0; i_r < grid_.numberSmootherCircles(); i_r++) {
            applyCircleSection(i_r, result, x);
        }
        for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) {
            applyRadialSection(i_theta, result, x);
        }
    }
    /* Multi-threaded execution */
    else {
        const int num_circle_tasks        = grid_.numberSmootherCircles();
        const int additional_radial_tasks = grid_.ntheta() % 3;
        const int num_radial_tasks        = grid_.ntheta() - additional_radial_tasks;

        #pragma omp parallel num_threads(num_omp_threads_)
        {
            /* Circle Section 0 */
            #pragma omp for
            for (int circle_task = 0; circle_task < num_circle_tasks; circle_task += 3) {
                int i_r = grid_.numberSmootherCircles() - circle_task - 1;
                applyCircleSection(i_r, result, x);
            }
            /* Circle Section 1 */
            #pragma omp for
            for (int circle_task = 1; circle_task < num_circle_tasks; circle_task += 3) {
                int i_r = grid_.numberSmootherCircles() - circle_task - 1;
                applyCircleSection(i_r, result, x);
            }
            /* Circle Section 2 */
            #pragma omp for nowait
            for (int circle_task = 2; circle_task < num_circle_tasks; circle_task += 3) {
                int i_r = grid_.numberSmootherCircles() - circle_task - 1;
                applyCircleSection(i_r, result, x);
            }

            /* Radial Section 0 */
            #pragma omp for
            for (int radial_task = 0; radial_task < num_radial_tasks; radial_task += 3) {
                if (radial_task > 0) {
                    int i_theta = radial_task + additional_radial_tasks;
                    applyRadialSection(i_theta, result, x);
                }
                else {
                    if (additional_radial_tasks == 0) {
                        applyRadialSection(0, result, x);
                    }
                    else if (additional_radial_tasks >= 1) {
                        applyRadialSection(0, result, x);
                        applyRadialSection(1, result, x);
                    }
                }
            }
            /* Radial Section 1 */
            #pragma omp for
            for (int radial_task = 1; radial_task < num_radial_tasks; radial_task += 3) {
                if (radial_task > 1) {
                    int i_theta = radial_task + additional_radial_tasks;
                    applyRadialSection(i_theta, result, x);
                }
                else {
                    if (additional_radial_tasks == 0) {
                        applyRadialSection(1, result, x);
                    }
                    else if (additional_radial_tasks == 1) {
                        applyRadialSection(2, result, x);
                    }
                    else if (additional_radial_tasks == 2) {
                        applyRadialSection(2, result, x);
                        applyRadialSection(3, result, x);
                    }
                }
            }
            /* Radial Section 2 */
            #pragma omp for
            for (int radial_task = 2; radial_task < num_radial_tasks; radial_task += 3) {
                int i_theta = radial_task + additional_radial_tasks;
                applyRadialSection(i_theta, result, x);
            }
        }
    }
}
// clang-format on