source: src/samples/discretization_poisson_fv.cpp@ dfed1c

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

Major vmg update.

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