source: src/samples/discretization_poisson_fv.cpp@ a72216

Last change on this file since a72216 was a72216, checked in by Olaf Lenz <olenz@…>, 13 years ago

Fixed permissions.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@2428 5161e1c8-67bf-11de-9fd5-51895aff932f
  • Property mode set to 100644
File size: 8.8 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 discretization_poisson_fv.cpp
21 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
22 * @date Mon Apr 18 13:03:47 2011
23 *
24 * @brief Finite volume discretization for the Poisson
25 * equation. Absolutely equivalent to the finite
26 * difference discretization unless you use
27 * hierarchically coarsened grids.
28 *
29 */
30
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35#include "samples/discretization_poisson_fv.hpp"
36
37using namespace VMG;
38
39void DiscretizationPoissonFV::SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const
40{
41 Index i_c, i_f;
42
43 const Vector h2_inv = 0.5 / sol_f.Extent().MeshWidth();
44
45 const Index b1_c = sol_c.Local().FinerBegin();
46 const Index b2_c = sol_c.Local().FinerEnd() - 1;
47 const Index b1_f = 0;
48 const Index b2_f = rhs_f.Local().SizeTotal() - 1;
49
50 const Index begin_f = sol_f.Local().Begin();
51 const Index end_f = sol_f.Local().End();
52 const Index begin_c = sol_c.Local().FinerBegin();
53
54 const vmg_float c_1_3 = 1.0 / 3.0;
55 const Vector c_2_3_sp = 2.0 / 3.0 * sol_f.Extent().MeshWidth();
56 const vmg_float c_4_3 = 4.0 / 3.0;
57
58 //
59 // X-direction
60 //
61 for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
62 for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z()) {
63 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = (sol_c(b1_c.X()-1,i_c.Y(),i_c.Z()) - sol_f(b1_f.X()+1,i_f.Y(),i_f.Z())) * h2_inv.X();
64 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = (sol_c(b2_c.X()+1,i_c.Y(),i_c.Z()) - sol_f(b2_f.X()-1,i_f.Y(),i_f.Z())) * h2_inv.X();
65 }
66
67 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
68 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2) {
69 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()) + rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()));
70 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()) + rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()));
71 }
72
73 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
74 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2) {
75 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b1_f.X(),i_f.Y(),i_f.Z()-1) + rhs_f(b1_f.X(),i_f.Y(),i_f.Z()+1));
76 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b2_f.X(),i_f.Y(),i_f.Z()-1) + rhs_f(b2_f.X(),i_f.Y(),i_f.Z()+1));
77 }
78
79 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
80 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2) {
81
82 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.25 * (rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()-1) +
83 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()-1) +
84 rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()+1) +
85 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()+1));
86
87 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.25 * (rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()-1) +
88 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()-1) +
89 rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()+1) +
90 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()+1));
91 }
92
93 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
94 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z()) {
95
96 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = sol_f(b1_f.X(),i_f.Y(),i_f.Z()) = c_2_3_sp.X() * rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) +
97 c_4_3 * sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()) -
98 c_1_3 * sol_f(b1_f.X()+2,i_f.Y(),i_f.Z());
99
100 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = sol_f(b2_f.X(),i_f.Y(),i_f.Z()) = c_2_3_sp.X() * rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) +
101 c_4_3 * sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()) -
102 c_1_3 * sol_f(b2_f.X()-2,i_f.Y(),i_f.Z());
103 }
104
105 //
106 // Y-direction
107 //
108
109 for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
110 for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z()) {
111 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = (sol_c(i_c.X(),b1_c.Y()-1,i_c.Z()) - sol_f(i_f.X(),b1_f.Y()+1,i_f.Z())) * h2_inv.Y();
112 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = (sol_c(i_c.X(),b2_c.Y()+1,i_c.Z()) - sol_f(i_f.X(),b2_f.Y()-1,i_f.Z())) * h2_inv.Y();
113 }
114
115 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
116 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2) {
117 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()));
118 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b2_f.Z(),i_f.Z()));
119 }
120
121 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
122 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2) {
123 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X(),b1_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b1_f.Y(),i_f.Z()+1));
124 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X(),b2_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b2_f.Y(),i_f.Z()+1));
125 }
126
127 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
128 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2) {
129
130 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()-1) +
131 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()-1) +
132 rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()+1) +
133 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()+1));
134
135 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()-1) +
136 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()-1) +
137 rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()+1) +
138 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()+1));
139 }
140
141 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
142 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z()) {
143
144 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = sol_f(i_f.X(),b1_f.Y(),i_f.Z()) = c_2_3_sp.Y() * rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) +
145 c_4_3 * sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()) -
146 c_1_3 * sol_f(i_f.X(),b1_f.Y()+2,i_f.Z());
147
148 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = sol_f(i_f.X(),b2_f.Y(),i_f.Z()) = c_2_3_sp.Y() * rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) +
149 c_4_3 * sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()) -
150 c_1_3 * sol_f(i_f.X(),b2_f.Y()-2,i_f.Z());
151 }
152
153 //
154 // Z-direction
155 //
156
157 for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
158 for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y()) {
159 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = (sol_c(i_c.X(),i_c.Y(),b1_c.Z()-1) - sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1)) * h2_inv.Z();
160 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = (sol_c(i_c.X(),i_c.Y(),b2_c.Z()+1) - sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1)) * h2_inv.Z();
161 }
162
163 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
164 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2) {
165 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,i_f.Y(),b1_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b1_f.Z()));
166 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,i_f.Y(),b2_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b2_f.Z()));
167 }
168
169 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
170 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2) {
171 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (rhs_f(i_f.X(),i_f.Y()-1,b1_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b1_f.Z()));
172 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (rhs_f(i_f.X(),i_f.Y()-1,b2_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b2_f.Z()));
173 }
174
175 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
176 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2) {
177 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b1_f.Z()) +
178 rhs_f(i_f.X()+1,i_f.Y()-1,b1_f.Z()) +
179 rhs_f(i_f.X()-1,i_f.Y()+1,b1_f.Z()) +
180 rhs_f(i_f.X()+1,i_f.Y()+1,b1_f.Z()));
181
182 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b2_f.Z()) +
183 rhs_f(i_f.X()+1,i_f.Y()-1,b2_f.Z()) +
184 rhs_f(i_f.X()-1,i_f.Y()+1,b2_f.Z()) +
185 rhs_f(i_f.X()+1,i_f.Y()+1,b2_f.Z()));
186 }
187
188 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
189 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y()) {
190
191 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = sol_f(i_f.X(),i_f.Y(),b1_f.Z()) = c_2_3_sp.Z() * rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) +
192 c_4_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1) -
193 c_1_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+2);
194
195 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = sol_f(i_f.X(),i_f.Y(),b2_f.Z()) = c_2_3_sp.Z() * rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) +
196 c_4_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1) -
197 c_1_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-2);
198 }
199
200#ifdef DEBUG_MATRIX_CHECKS
201 rhs_f.IsConsistent();
202 sol_f.IsConsistent();
203#endif
204}
Note: See TracBrowser for help on using the repository browser.