source: ThirdParty/vmg/src/smoother/gsrb_poisson_2.cpp@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 7faa5c, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit 'de061d9d851257a04e924d4472df4523d33bb08b' as 'ThirdParty/vmg'
  • Property mode set to 100644
File size: 4.9 KB
Line 
1/*
2 * vmg - a versatile multigrid solver
3 * Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
4 *
5 * vmg is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * vmg is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/**
20 * @file gsrb_poisson_2.cpp
21 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
22 * @date Fri May 11 18:30:20 2012
23 *
24 * @brief Gauss-Seidel Red Black method, specialized to
25 * the Poisson equation. Performance improved by
26 * explicit loop unrolling.
27 *
28 */
29
30#ifdef HAVE_CONFIG_H
31#include <libvmg_config.h>
32#endif
33
34#ifdef HAVE_MPI
35#include <mpi.h>
36#endif
37
38#include "base/helper.hpp"
39#include "comm/comm.hpp"
40#include "grid/grid.hpp"
41#include "smoother/gsrb_poisson_2.hpp"
42#include "mg.hpp"
43
44using namespace VMG;
45
46static inline void ComputePartial(Grid& sol, Grid& rhs,
47 const Index& begin, const Index& end,
48 const vmg_float& prefactor, const int& off)
49{
50 const vmg_float fac = 1.0 / 6.0;
51
52 for (int i=begin.X(); i<end.X(); ++i)
53 for (int j=begin.Y(); j<end.Y(); ++j)
54 for (int k=begin.Z() + (i + j + begin.Z() + off) % 2; k<end.Z(); k+=2)
55 sol(i,j,k) = prefactor * rhs.GetVal(i,j,k) + fac * (sol.GetVal(i-1,j ,k ) +
56 sol.GetVal(i+1,j ,k ) +
57 sol.GetVal(i ,j-1,k ) +
58 sol.GetVal(i ,j+1,k ) +
59 sol.GetVal(i ,j ,k-1) +
60 sol.GetVal(i ,j ,k+1));
61}
62
63void GaussSeidelRBPoisson2::Compute(Grid& sol, Grid& rhs)
64{
65 const vmg_float prefactor_inv = Helper::pow_2(sol.Extent().MeshWidth().Max()) / 6.0;
66 const int off = rhs.Global().LocalBegin().Sum() - rhs.Global().GlobalBegin().Sum() + rhs.Local().HaloSize1().Sum();
67 const LocalIndices& local = rhs.Local();
68 Comm& comm = *MG::GetComm();
69
70 /*
71 * Compute first halfstep
72 */
73
74 // Start asynchronous communication
75 comm.CommToGhostsAsyncStart(sol);
76
77 // Smooth part not depending on ghost cells
78 ComputePartial(sol, rhs,
79 local.Begin()+1, local.End()-1,
80 prefactor_inv, off+1);
81
82 // Finish asynchronous communication
83 comm.CommToGhostsAsyncFinish(sol);
84
85 /*
86 * Smooth near boundary cells
87 */
88
89 ComputePartial(sol, rhs,
90 local.Begin(),
91 Index(local.Begin().X()+1, local.End().Y(), local.End().Z()),
92 prefactor_inv, off+1);
93
94 ComputePartial(sol, rhs,
95 Index(local.End().X()-1, local.Begin().Y(), local.Begin().Z()),
96 local.End(),
97 prefactor_inv, off+1);
98
99 ComputePartial(sol, rhs,
100 Index(local.Begin().X()+1, local.Begin().Y(), local.Begin().Z()),
101 Index(local.End().X()-1, local.Begin().Y()+1, local.End().Z()),
102 prefactor_inv, off+1);
103
104 ComputePartial(sol, rhs,
105 Index(local.Begin().X()+1, local.End().Y()-1, local.Begin().Z()),
106 Index(local.End().X()-1, local.End().Y(), local.End().Z()),
107 prefactor_inv, off+1);
108
109 ComputePartial(sol, rhs,
110 Index(local.Begin().X()+1, local.Begin().Y()+1, local.Begin().Z()),
111 Index(local.End().X()-1, local.End().Y()-1, local.Begin().Z()+1),
112 prefactor_inv, off+1);
113
114 ComputePartial(sol, rhs,
115 Index(local.Begin().X()+1, local.Begin().Y()+1, local.End().Z()-1),
116 Index(local.End().X()-1, local.End().Y()-1, local.End().Z()),
117 prefactor_inv, off+1);
118
119 /*
120 * Compute second halfstep
121 */
122
123 // Start asynchronous communication
124 comm.CommToGhostsAsyncStart(sol);
125
126 // Smooth part not depending on ghost cells
127 ComputePartial(sol, rhs,
128 local.Begin()+1, local.End()-1,
129 prefactor_inv, off);
130
131 // Finish asynchronous communication
132 comm.CommToGhostsAsyncFinish(sol);
133
134 /*
135 * Smooth near boundary cells
136 */
137
138 ComputePartial(sol, rhs,
139 local.Begin(),
140 Index(local.Begin().X()+1, local.End().Y(), local.End().Z()),
141 prefactor_inv, off);
142
143 ComputePartial(sol, rhs,
144 Index(local.End().X()-1, local.Begin().Y(), local.Begin().Z()),
145 local.End(),
146 prefactor_inv, off);
147
148 ComputePartial(sol, rhs,
149 Index(local.Begin().X()+1, local.Begin().Y(), local.Begin().Z()),
150 Index(local.End().X()-1, local.Begin().Y()+1, local.End().Z()),
151 prefactor_inv, off);
152
153 ComputePartial(sol, rhs,
154 Index(local.Begin().X()+1, local.End().Y()-1, local.Begin().Z()),
155 Index(local.End().X()-1, local.End().Y(), local.End().Z()),
156 prefactor_inv, off);
157
158 ComputePartial(sol, rhs,
159 Index(local.Begin().X()+1, local.Begin().Y()+1, local.Begin().Z()),
160 Index(local.End().X()-1, local.End().Y()-1, local.Begin().Z()+1),
161 prefactor_inv, off);
162
163 ComputePartial(sol, rhs,
164 Index(local.Begin().X()+1, local.Begin().Y()+1, local.End().Z()-1),
165 Index(local.End().X()-1, local.End().Y()-1, local.End().Z()),
166 prefactor_inv, off);
167}
Note: See TracBrowser for help on using the repository browser.