source: src/samples/discretization_poisson_fv.cpp@ bd61b5c

Last change on this file since bd61b5c was 8180d8, checked in by Julian Iseringhausen <julian.iseringhausen@…>, 13 years ago

Merge stashed open boundary stuff.
  • Property mode set to 100644
File size: 11.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 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 "grid/grid_index_translations.hpp"
36#include "samples/discretization_poisson_fv.hpp"
37
38using namespace VMG;
39
40void DiscretizationPoissonFV::InitDiscretizationPoissonFV()
41{
42 switch (order)
43 {
44 case 2:
45 stencil.SetDiag(6.0);
46 stencil.push_back(-1, 0, 0, -1.0);
47 stencil.push_back( 1, 0, 0, -1.0);
48 stencil.push_back( 0, -1, 0, -1.0);
49 stencil.push_back( 0, 1, 0, -1.0);
50 stencil.push_back( 0, 0, -1, -1.0);
51 stencil.push_back( 0, 0, 1, -1.0);
52 break;
53 case 4:
54 stencil.SetDiag(24.0/6.0);
55 stencil.push_back(-1, 0, 0, -2.0/6.0);
56 stencil.push_back( 1, 0, 0, -2.0/6.0);
57 stencil.push_back( 0, -1, 0, -2.0/6.0);
58 stencil.push_back( 0, 1, 0, -2.0/6.0);
59 stencil.push_back( 0, 0, -1, -2.0/6.0);
60 stencil.push_back( 0, 0, 1, -2.0/6.0);
61 stencil.push_back(-1, -1, 0, -1.0/6.0);
62 stencil.push_back(-1, 1, 0, -1.0/6.0);
63 stencil.push_back( 1, -1, 0, -1.0/6.0);
64 stencil.push_back( 1, 1, 0, -1.0/6.0);
65 stencil.push_back(-1, 0, -1, -1.0/6.0);
66 stencil.push_back(-1, 0, 1, -1.0/6.0);
67 stencil.push_back( 1, 0, -1, -1.0/6.0);
68 stencil.push_back( 1, 0, 1, -1.0/6.0);
69 stencil.push_back( 0, -1, -1, -1.0/6.0);
70 stencil.push_back( 0, -1, 1, -1.0/6.0);
71 stencil.push_back( 0, 1, -1, -1.0/6.0);
72 stencil.push_back( 0, 1, 1, -1.0/6.0);
73 break;
74 default:
75 assert(0 != "vmg choose discretization order 2 or 4");
76 break;
77 }
78}
79
80void DiscretizationPoissonFV::ModifyRightHandSide()
81{
82 if (order == 4) {
83
84 Grid& rhs = MG::GetRhsMaxLevel();
85
86 Stencil stencil(6.0/12.0);
87 stencil.push_back(-1, 0, 0, 1.0/12.0);
88 stencil.push_back( 1, 0, 0, 1.0/12.0);
89 stencil.push_back( 0, -1, 0, 1.0/12.0);
90 stencil.push_back( 0, 1, 0, 1.0/12.0);
91 stencil.push_back( 0, 0, -1, 1.0/12.0);
92 stencil.push_back( 0, 0, 1, 1.0/12.0);
93
94 stencil.Apply(rhs);
95
96 }
97}
98
99void DiscretizationPoissonFV::SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const
100{
101/*
102 //FIXME
103 Index i_c, i_f;
104
105 Index b1_c, b2_c, b1_f, b2_f;
106 GridIndexTranslations::GetInnerBoundary(sol_c, b1_c, b2_c, sol_f, b1_f, b2_f);
107*/
108/*
109 const Index b1_c = sol_c.Local().FinerBegin();
110 const Index b2_c = sol_c.Local().FinerEnd() - 1;
111 const Index b1_f = 0;
112 const Index b2_f = rhs_f.Local().SizeTotal() - 1;
113*/
114
115/*
116 const Index begin_f = sol_f.Local().Begin();
117 const Index end_f = sol_f.Local().End();
118 const Index begin_c = sol_c.Local().FinerBegin();
119
120 const Vector h2_inv = 0.5 / sol_f.Extent().MeshWidth();
121 const vmg_float c_1_3 = 1.0 / 3.0;
122 const Vector c_2_3_sp = 2.0 / 3.0 * sol_f.Extent().MeshWidth();
123 const vmg_float c_4_3 = 4.0 / 3.0;
124
125 Helper::AssertVectorsEqual(sol_c.GetSpatialPos(sol_c.Local().FinerBegin()), sol_f.GetSpatialPos(sol_f.Local().Begin()));
126 Helper::AssertVectorsEqual(sol_c.GetSpatialPos(sol_c.Local().FinerEnd()-1), sol_f.GetSpatialPos(sol_f.Local().End()-1));
127
128 //
129 // X-direction
130 //
131 i_f = begin_f; i_c = begin_c;
132 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())
133 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()) {
134 Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
135 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();
136 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();
137 }
138
139 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
140 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2) {
141 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()));
142 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()));
143 }
144
145 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
146 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2) {
147 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));
148 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));
149 }
150
151 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
152 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2) {
153
154 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) +
155 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()-1) +
156 rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()+1) +
157 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()+1));
158
159 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) +
160 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()-1) +
161 rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()+1) +
162 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()+1));
163 }
164
165 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
166 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z()) {
167
168 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()) +
169 c_4_3 * sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()) -
170 c_1_3 * sol_f(b1_f.X()+2,i_f.Y(),i_f.Z());
171
172 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()) +
173 c_4_3 * sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()) -
174 c_1_3 * sol_f(b2_f.X()-2,i_f.Y(),i_f.Z());
175 }
176
177 //
178 // Y-direction
179 //
180 i_f = begin_f; i_c = begin_c;
181 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())
182 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()) {
183 Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
184 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();
185 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();
186 }
187
188 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
189 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2) {
190 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()));
191 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()));
192 }
193
194 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
195 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2) {
196 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));
197 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));
198 }
199
200 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
201 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2) {
202
203 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) +
204 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()-1) +
205 rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()+1) +
206 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()+1));
207
208 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) +
209 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()-1) +
210 rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()+1) +
211 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()+1));
212 }
213
214 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
215 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z()) {
216
217 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()) +
218 c_4_3 * sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()) -
219 c_1_3 * sol_f(i_f.X(),b1_f.Y()+2,i_f.Z());
220
221 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()) +
222 c_4_3 * sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()) -
223 c_1_3 * sol_f(i_f.X(),b2_f.Y()-2,i_f.Z());
224 }
225
226 //
227 // Z-direction
228 //
229 i_f = begin_f; i_c = begin_c;
230 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())
231 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()) {
232 Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
233 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();
234 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();
235 }
236
237 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
238 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2) {
239 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()));
240 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()));
241 }
242
243 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
244 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2) {
245 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()));
246 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()));
247 }
248
249 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
250 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2) {
251 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()) +
252 rhs_f(i_f.X()+1,i_f.Y()-1,b1_f.Z()) +
253 rhs_f(i_f.X()-1,i_f.Y()+1,b1_f.Z()) +
254 rhs_f(i_f.X()+1,i_f.Y()+1,b1_f.Z()));
255
256 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()) +
257 rhs_f(i_f.X()+1,i_f.Y()-1,b2_f.Z()) +
258 rhs_f(i_f.X()-1,i_f.Y()+1,b2_f.Z()) +
259 rhs_f(i_f.X()+1,i_f.Y()+1,b2_f.Z()));
260 }
261
262 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
263 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y()) {
264
265 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()) +
266 c_4_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1) -
267 c_1_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+2);
268
269 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()) +
270 c_4_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1) -
271 c_1_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-2);
272 }
273*/
274#ifdef DEBUG_MATRIX_CHECKS
275 rhs_f.IsConsistent();
276 sol_f.IsConsistent();
277#endif
278}
Note: See TracBrowser for help on using the repository browser.