source: src/samples/discretization_poisson_fv.cpp@ ac6d04

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

Merge recent changes of the vmg library into ScaFaCos.

Includes a fix for the communication problems on Jugene.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1666 5161e1c8-67bf-11de-9fd5-51895aff932f
  • Property mode set to 100644
File size: 8.0 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 Vector h2_inv = 0.5 / sol_f.Extent().MeshWidth();
26
27 const Index b1_c = sol_c.Local().FinerBeginFoo();
28 const Index b2_c = sol_c.Local().FinerEndFoo() - 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().FinerBeginFoo();
35
36 const vmg_float c_1_3 = 1.0 / 3.0;
37 const Vector c_2_3_sp = 2.0 / 3.0 * sol_f.Extent().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 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();
46 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();
47 }
48
49 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
50 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2) {
51 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()));
52 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()));
53 }
54
55 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
56 for (i_f.Z()=begin_f.Z()+1; 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(),i_f.Z()-1) + rhs_f(b1_f.X(),i_f.Y(),i_f.Z()+1));
58 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));
59 }
60
61 for (i_f.Y()=begin_f.Y()+1; 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
64 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) +
65 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()-1) +
66 rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()+1) +
67 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()+1));
68
69 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) +
70 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()-1) +
71 rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()+1) +
72 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()+1));
73 }
74
75 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
76 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z()) {
77
78 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()) +
79 c_4_3 * sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()) -
80 c_1_3 * sol_f(b1_f.X()+2,i_f.Y(),i_f.Z());
81
82 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()) +
83 c_4_3 * sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()) -
84 c_1_3 * sol_f(b2_f.X()-2,i_f.Y(),i_f.Z());
85 }
86
87 //
88 // Y-direction
89 //
90
91 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())
92 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()) {
93 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();
94 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();
95 }
96
97 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
98 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2) {
99 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()));
100 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()));
101 }
102
103 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
104 for (i_f.Z()=begin_f.Z()+1; 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(),b1_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b1_f.Y(),i_f.Z()+1));
106 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));
107 }
108
109 for (i_f.X()=begin_f.X()+1; 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
112 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) +
113 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()-1) +
114 rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()+1) +
115 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()+1));
116
117 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) +
118 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()-1) +
119 rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()+1) +
120 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()+1));
121 }
122
123 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
124 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z()) {
125
126 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()) +
127 c_4_3 * sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()) -
128 c_1_3 * sol_f(i_f.X(),b1_f.Y()+2,i_f.Z());
129
130 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()) +
131 c_4_3 * sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()) -
132 c_1_3 * sol_f(i_f.X(),b2_f.Y()-2,i_f.Z());
133 }
134
135 //
136 // Z-direction
137 //
138
139 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())
140 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()) {
141 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();
142 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();
143 }
144
145 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
146 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2) {
147 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()));
148 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()));
149 }
150
151 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
152 for (i_f.Y()=begin_f.Y()+1; 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(),i_f.Y()-1,b1_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b1_f.Z()));
154 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()));
155 }
156
157 for (i_f.X()=begin_f.X()+1; 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.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b1_f.Z()) +
160 rhs_f(i_f.X()+1,i_f.Y()-1,b1_f.Z()) +
161 rhs_f(i_f.X()-1,i_f.Y()+1,b1_f.Z()) +
162 rhs_f(i_f.X()+1,i_f.Y()+1,b1_f.Z()));
163
164 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()) +
165 rhs_f(i_f.X()+1,i_f.Y()-1,b2_f.Z()) +
166 rhs_f(i_f.X()-1,i_f.Y()+1,b2_f.Z()) +
167 rhs_f(i_f.X()+1,i_f.Y()+1,b2_f.Z()));
168 }
169
170 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
171 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y()) {
172
173 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()) +
174 c_4_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1) -
175 c_1_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+2);
176
177 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()) +
178 c_4_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1) -
179 c_1_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-2);
180 }
181
182#ifdef DEBUG_MATRIX_CHECKS
183 rhs_f.IsConsistent();
184 sol_f.IsConsistent();
185#endif
186}
Note: See TracBrowser for help on using the repository browser.