source: ThirdParty/vmg/src/smoother/gsrb.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: 5.5 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.cpp
21 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
22 * @date Mon Apr 18 13:08:20 2011
23 *
24 * @brief Gauss-Seidel Red Black method
25 *
26 */
27
28#ifdef HAVE_CONFIG_H
29#include <libvmg_config.h>
30#endif
31
32#ifdef HAVE_MPI
33#include <mpi.h>
34#endif
35
36#include "base/discretization.hpp"
37#include "base/stencil.hpp"
38#include "comm/comm.hpp"
39#include "grid/grid.hpp"
40#include "smoother/gsrb.hpp"
41
42using namespace VMG;
43
44static inline void ComputePartial(Grid& sol, Grid& rhs, const Stencil& mat,
45 const Index& begin, const Index& end,
46 const vmg_float& prefactor, const vmg_float& diag_inv,
47 const int& off)
48{
49 int i,j,k;
50 vmg_float temp;
51 Stencil::iterator iter;
52
53 for (i=begin.X(); i<end.X(); ++i)
54 for (j=begin.Y(); j<end.Y(); ++j) {
55
56 int z_begin = begin.Z() + (i + j + begin.Z() + off) % 2;
57
58#ifdef DEBUG
59 int off_sum = MG::GetComm()->LevelSum(rhs, z_begin - begin.Z());
60 assert(z_begin - begin.Z() == 0 || z_begin - begin.Z() == 1);
61 assert(off_sum == 0 || off_sum == MG::GetComm()->Size(rhs));
62#endif /* DEBUG */
63
64 for (k=z_begin; k<end.Z(); k+=2) {
65
66 temp = prefactor * rhs.GetVal(i,j,k);
67
68 for (iter=mat.begin(); iter!=mat.end(); ++iter)
69 temp -= iter->Val() * sol.GetVal(i+iter->Disp().X(),
70 j+iter->Disp().Y(),
71 k+iter->Disp().Z());
72
73 sol(i,j,k) = temp * diag_inv;
74
75 }
76 }
77}
78
79void GaussSeidelRB::Compute(Grid& sol, Grid& rhs)
80{
81 const Stencil& mat = MG::GetDiscretization()->GetStencil();
82 const vmg_float prefactor_inv = 1.0 / MG::GetDiscretization()->OperatorPrefactor(sol);
83 const vmg_float diag_inv = 1.0 / mat.GetDiag();
84 const int off = rhs.Global().LocalBegin().Sum() - rhs.Global().GlobalBegin().Sum() + rhs.Local().HaloSize1().Sum();
85 const LocalIndices& local = rhs.Local();
86 Comm& comm = *MG::GetComm();
87
88 /*
89 * Compute first halfstep
90 */
91
92 // Start asynchronous communication
93 comm.CommToGhostsAsyncStart(sol);
94
95 // Smooth part not depending on ghost cells
96 ComputePartial(sol, rhs, mat,
97 local.Begin()+1, local.End()-1,
98 prefactor_inv, diag_inv, off+1);
99
100 // Finish asynchronous communication
101 comm.CommToGhostsAsyncFinish(sol);
102
103 /*
104 * Smooth near boundary cells
105 */
106
107 ComputePartial(sol, rhs, mat,
108 local.Begin(),
109 Index(local.Begin().X()+1, local.End().Y(), local.End().Z()),
110 prefactor_inv, diag_inv, off+1);
111
112 ComputePartial(sol, rhs, mat,
113 Index(local.End().X()-1, local.Begin().Y(), local.Begin().Z()),
114 local.End(),
115 prefactor_inv, diag_inv, off+1);
116
117 ComputePartial(sol, rhs, mat,
118 Index(local.Begin().X()+1, local.Begin().Y(), local.Begin().Z()),
119 Index(local.End().X()-1, local.Begin().Y()+1, local.End().Z()),
120 prefactor_inv, diag_inv, off+1);
121
122 ComputePartial(sol, rhs, mat,
123 Index(local.Begin().X()+1, local.End().Y()-1, local.Begin().Z()),
124 Index(local.End().X()-1, local.End().Y(), local.End().Z()),
125 prefactor_inv, diag_inv, off+1);
126
127 ComputePartial(sol, rhs, mat,
128 Index(local.Begin().X()+1, local.Begin().Y()+1, local.Begin().Z()),
129 Index(local.End().X()-1, local.End().Y()-1, local.Begin().Z()+1),
130 prefactor_inv, diag_inv, off+1);
131
132 ComputePartial(sol, rhs, mat,
133 Index(local.Begin().X()+1, local.Begin().Y()+1, local.End().Z()-1),
134 Index(local.End().X()-1, local.End().Y()-1, local.End().Z()),
135 prefactor_inv, diag_inv, off+1);
136
137 /*
138 * Compute second halfstep
139 */
140
141 // Start asynchronous communication
142 comm.CommToGhostsAsyncStart(sol);
143
144 // Smooth part not depending on ghost cells
145 ComputePartial(sol, rhs, mat,
146 local.Begin()+1, local.End()-1,
147 prefactor_inv, diag_inv, off);
148
149 // Finish asynchronous communication
150 comm.CommToGhostsAsyncFinish(sol);
151
152 /*
153 * Smooth near boundary cells
154 */
155
156 ComputePartial(sol, rhs, mat,
157 local.Begin(),
158 Index(local.Begin().X()+1, local.End().Y(), local.End().Z()),
159 prefactor_inv, diag_inv, off);
160
161 ComputePartial(sol, rhs, mat,
162 Index(local.End().X()-1, local.Begin().Y(), local.Begin().Z()),
163 local.End(),
164 prefactor_inv, diag_inv, off);
165
166 ComputePartial(sol, rhs, mat,
167 Index(local.Begin().X()+1, local.Begin().Y(), local.Begin().Z()),
168 Index(local.End().X()-1, local.Begin().Y()+1, local.End().Z()),
169 prefactor_inv, diag_inv, off);
170
171 ComputePartial(sol, rhs, mat,
172 Index(local.Begin().X()+1, local.End().Y()-1, local.Begin().Z()),
173 Index(local.End().X()-1, local.End().Y(), local.End().Z()),
174 prefactor_inv, diag_inv, off);
175
176 ComputePartial(sol, rhs, mat,
177 Index(local.Begin().X()+1, local.Begin().Y()+1, local.Begin().Z()),
178 Index(local.End().X()-1, local.End().Y()-1, local.Begin().Z()+1),
179 prefactor_inv, diag_inv, off);
180
181 ComputePartial(sol, rhs, mat,
182 Index(local.Begin().X()+1, local.Begin().Y()+1, local.End().Z()-1),
183 Index(local.End().X()-1, local.End().Y()-1, local.End().Z()),
184 prefactor_inv, diag_inv, off);
185}
Note: See TracBrowser for help on using the repository browser.