source: src/samples/discretization_poisson_fv.cpp@ f57182

Last change on this file since f57182 was f57182, checked in by Julian Iseringhausen <isering@…>, 13 years ago

Open boundary conditions.

Conflicts:

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