AMR-Solver  1.0
Block-based Octree AMR grid flow solver
poisson.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <string>
3 #include "octreegrid.h"
4 #include "poisson.h"
5 #include "ghost.h"
6 #include <iostream>
7 #include <cmath>
8 
9 using myOctree::Field;
10 using myOctree::VecField;
11 using myOctree::Octree;
12 
13 using myOctree::nx_block;
14 using myOctree::ny_block;
15 using myOctree::nz_block;
16 using myOctree::pad;
17 
18 using myOctree::level_nodes;
19 
20 namespace amrsolver {
21 
22 //
23 void jacobi(int level, std::string name) {
24 
25  myOctree::create_lists_of_level_nodes();
26 
27  if(myOctree::level_nodes[level].empty()) {
28  std::cerr << "No blocks in level " << level << std::endl;
29  exit(1);
30  }
31 
32  double global_res = 1.0;
33 
34  long int loop = 0;
35 
36  //iteration loop
37  while(loop <250 ) {
38  //while(global_res > pow(10,-7)) {
39 
40  global_res = 0.0;
41 
42  //exchange ghost values
43  myOctree::exchange_ghost_val(level, name);
44 
45  //loop over all the nodes
46  for (std::list<Octree*>::iterator it = level_nodes[level].begin(), end = level_nodes[level].end(); it != end; ++it) {
47 
48  for(int l = 0; l<myOctree::scalar_fields.size() ; l++) {
49 
50  Field *f = (*it)->get_block_data()->scalarfields[l];
51  // std::cerr << f->name << std::endl;
52 
53  if( f->name == name ) {
54 
55  global_res = jacobi_for_field(*it, f, global_res);
56 
57  }
58  }
59  }
60 
61  std::cout << "Residual " << global_res << std::endl;
62 
63  loop++;
64 
65  }
66 }
67 
69 void gauss_seidel(int level, std::string name) {
70 
71  myOctree::create_lists_of_level_nodes();
72 
73  if(myOctree::level_nodes[level].empty()) {
74  std::cerr << "No blocks in level " << level << std::endl;
75  exit(1);
76  }
77 
78  double global_res = 1.0;
79 
80  long int loop = 0;
81 
82  //iteration loop
83  //while(loop <250 ) {
84  while(global_res > pow(10,-7)) {
85 
86  global_res = 0.0;
87 
88  //exchange ghost values
89  myOctree::exchange_ghost_val(level, name);
90 
91  //loop over all the nodes
92  for (std::list<Octree*>::iterator it = level_nodes[level].begin(), end = level_nodes[level].end(); it != end; ++it) {
93 
94  for(int l = 0; l<myOctree::scalar_fields.size() ; l++) {
95 
96  Field *f = (*it)->get_block_data()->scalarfields[l];
97  if( f->name == name ) {
98 
99  global_res = gauss_seidel_red(*it, f, global_res);
100 
101  }
102  }
103  }
104 
105  //exchange ghost values
106  myOctree::exchange_ghost_val(level, name);
107 
108  //loop over all the nodes
109  for (std::list<Octree*>::iterator it = level_nodes[level].begin(), end = level_nodes[level].end(); it != end; ++it) {
110 
111  for(int l = 0; l<myOctree::scalar_fields.size() ; l++) {
112 
113  Field *f = (*it)->get_block_data()->scalarfields[l];
114  if( f->name == name ) {
115 
116  global_res = gauss_seidel_black(*it, f, global_res);
117 
118  }
119  }
120  }
121 
122  std::cout << "Residual " << global_res << std::endl;
123 
124  loop++;
125 
126  }
127  std::cout << "No of loops " << loop << std::endl;
128 }
129 
130 double jacobi_for_field(Octree* node, Field* f, double global_res) {
131 
132  double res, err;
133  double dx, dy, dz;
134  double ue, uw, un, us, ut, ub, up;
135  double force = 1.0;
136  res = 0.0;
137 
138  Field temp(nx_block+2*pad,ny_block+2*pad,nz_block+2*pad, "temp");
139  //temp.set_field(0.0);
140 
141  dx = node->get_block_data()->dx;
142  dy = node->get_block_data()->dy;
143  dz = node->get_block_data()->dz;
144 
145 
146  //jocobi work - add here
147  for(int i=pad; i<nx_block+pad; i++) {
148  for(int j=pad; j<ny_block+pad; j++) {
149  for(int k=pad; k<nz_block+pad; k++) {
150 
151  ue = f->val[i+1][j][k];
152  uw = f->val[i-1][j][k];
153  un = f->val[i][j+1][k];
154  us = f->val[i][j-1][k];
155  ut = f->val[i][j][k+1];
156  ub = f->val[i][j][k-1];
157  up = f->val[i][j][k];
158 
159 
160  if(i==pad) temp.val[i][j][k] += (2.0*uw - up)/pow(dx,2);
161  else temp.val[i][j][k] += uw/pow(dx,2);
162  if(i==nx_block+pad-1) temp.val[i][j][k] += (2.0*ue - up)/pow(dx,2);
163  else temp.val[i][j][k] += ue/pow(dx,2);
164  if(j==pad) temp.val[i][j][k] += (2.0*us - up)/pow(dy,2);
165  else temp.val[i][j][k] += us/pow(dy,2);
166  if(j==ny_block+pad-1) temp.val[i][j][k] += (2.0*un - up)/pow(dy,2);
167  else temp.val[i][j][k] += un/pow(dy,2);
168  if(k==pad) temp.val[i][j][k] += (2.0*ub - up)/pow(dz,2);
169  else temp.val[i][j][k] += ub/pow(dz,2);
170  if(k==nz_block+pad-1) temp.val[i][j][k] += (2.0*ut - up)/pow(dz,2);
171  else temp.val[i][j][k] += ut/pow(dz,2);
172 
173  temp.val[i][j][k] -= force;
174  temp.val[i][j][k] /= 2.0;
175  temp.val[i][j][k] /= (1.0/pow(dx,2) + 1.0/pow(dy,2) + 1.0/pow(dz,2));
176 
177  // std::cerr << temp.val[i][j][k] << std::endl;
178 
179  /*magnitude to be considered*/
180  err = std::abs(temp.val[i][j][k] - f->val[i][j][k]);
181  if(err>res) res = err;
182 
183 
184  }
185  }
186  }
187 
188  //Transferring values from temp to field
189  for(int i=pad; i<nx_block+pad; i++) {
190  for(int j=pad; j<ny_block+pad; j++) {
191  for(int k=pad; k<nz_block+pad; k++) {
192 
193  f->val[i][j][k] = temp.val[i][j][k];
194 
195  }
196  }
197  }
198 
199 
200  if(res>global_res) global_res = res;
201 
202  return global_res;
203 
204 }
205 
207 double gauss_seidel_red(Octree* node, Field* f, double global_res) {
208 
209  double res, err;
210  double dx, dy, dz;
211  double ue, uw, un, us, ut, ub, up;
212  double prev_val, curr_val;
213  double force = 100.0;
214  double relax = 1.5;
215  res = 0.0;
216 
217  Field temp(nx_block+2*pad,ny_block+2*pad,nz_block+2*pad, "temp");
218  //temp.set_field(0.0);
219 
220  dx = node->get_block_data()->dx;
221  dy = node->get_block_data()->dy;
222  dz = node->get_block_data()->dz;
223 
224 
225  //jocobi work - add here
226  for(int i=pad; i<nx_block+pad; i++) {
227  for(int j=pad; j<ny_block+pad; j++) {
228  for(int k=pad; k<nz_block+pad; k++) {
229 
230 
231  if((i+j+k)%2==0) {
232 
233  ue = f->val[i+1][j][k];
234  uw = f->val[i-1][j][k];
235  un = f->val[i][j+1][k];
236  us = f->val[i][j-1][k];
237  ut = f->val[i][j][k+1];
238  ub = f->val[i][j][k-1];
239  up = f->val[i][j][k];
240 
241 
242  if(i==pad) temp.val[i][j][k] += (2.0*uw - up)/pow(dx,2);
243  else temp.val[i][j][k] += uw/pow(dx,2);
244  if(i==nx_block+pad-1) temp.val[i][j][k] += (2.0*ue - up)/pow(dx,2);
245  else temp.val[i][j][k] += ue/pow(dx,2);
246  if(j==pad) temp.val[i][j][k] += (2.0*us - up)/pow(dy,2);
247  else temp.val[i][j][k] += us/pow(dy,2);
248  if(j==ny_block+pad-1) temp.val[i][j][k] += (2.0*un - up)/pow(dy,2);
249  else temp.val[i][j][k] += un/pow(dy,2);
250  if(k==pad) temp.val[i][j][k] += (2.0*ub - up)/pow(dz,2);
251  else temp.val[i][j][k] += ub/pow(dz,2);
252  if(k==nz_block+pad-1) temp.val[i][j][k] += (2.0*ut - up)/pow(dz,2);
253  else temp.val[i][j][k] += ut/pow(dz,2);
254 
255  temp.val[i][j][k] -= force;
256  temp.val[i][j][k] /= 2.0;
257  temp.val[i][j][k] /= (1.0/pow(dx,2) + 1.0/pow(dy,2) + 1.0/pow(dz,2));
258 
259  prev_val = f->val[i][j][k];
260 
261  f->val[i][j][k] = (1.0 - relax) * f->val[i][j][k] + relax * temp.val[i][j][k];
262 
263  curr_val = f->val[i][j][k];
264 
265  /*magnitude to be considered*/
266  err = std::abs(curr_val - prev_val);
267  if(err>res) res = err;
268 
269  }
270 
271  }
272  }
273  }
274 
275 
276  if(res>global_res) global_res = res;
277 
278  return global_res;
279 
280 }
281 
283 double gauss_seidel_black(Octree* node, Field* f, double global_res) {
284 
285  double res, err;
286  double dx, dy, dz;
287  double ue, uw, un, us, ut, ub, up;
288  double prev_val, curr_val;
289  double force = 100.0;
290  double relax = 1.5;
291  res = 0.0;
292 
293  Field temp(nx_block+2*pad,ny_block+2*pad,nz_block+2*pad, "temp");
294  //temp.set_field(0.0);
295 
296  dx = node->get_block_data()->dx;
297  dy = node->get_block_data()->dy;
298  dz = node->get_block_data()->dz;
299 
300 
301  //jocobi work - add here
302  for(int i=pad; i<nx_block+pad; i++) {
303  for(int j=pad; j<ny_block+pad; j++) {
304  for(int k=pad; k<nz_block+pad; k++) {
305 
306  if((i+j+k)%2!=0) {
307 
308  ue = f->val[i+1][j][k];
309  uw = f->val[i-1][j][k];
310  un = f->val[i][j+1][k];
311  us = f->val[i][j-1][k];
312  ut = f->val[i][j][k+1];
313  ub = f->val[i][j][k-1];
314  up = f->val[i][j][k];
315 
316 
317  if(i==pad) temp.val[i][j][k] += (2.0*uw - up)/pow(dx,2);
318  else temp.val[i][j][k] += uw/pow(dx,2);
319  if(i==nx_block+pad-1) temp.val[i][j][k] += (2.0*ue - up)/pow(dx,2);
320  else temp.val[i][j][k] += ue/pow(dx,2);
321  if(j==pad) temp.val[i][j][k] += (2.0*us - up)/pow(dy,2);
322  else temp.val[i][j][k] += us/pow(dy,2);
323  if(j==ny_block+pad-1) temp.val[i][j][k] += (2.0*un - up)/pow(dy,2);
324  else temp.val[i][j][k] += un/pow(dy,2);
325  if(k==pad) temp.val[i][j][k] += (2.0*ub - up)/pow(dz,2);
326  else temp.val[i][j][k] += ub/pow(dz,2);
327  if(k==nz_block+pad-1) temp.val[i][j][k] += (2.0*ut - up)/pow(dz,2);
328  else temp.val[i][j][k] += ut/pow(dz,2);
329 
330  temp.val[i][j][k] -= force;
331  temp.val[i][j][k] /= 2.0;
332  temp.val[i][j][k] /= (1.0/pow(dx,2) + 1.0/pow(dy,2) + 1.0/pow(dz,2));
333 
334  prev_val = f->val[i][j][k];
335 
336  f->val[i][j][k] = (1.0 - relax) * f->val[i][j][k] + relax * temp.val[i][j][k];
337 
338  curr_val = f->val[i][j][k];
339 
340  /*magnitude to be considered*/
341  err = std::abs(curr_val - prev_val);
342  if(err>res) res = err;
343 
344  }
345  }
346  }
347  }
348 
349 
350  if(res>global_res) global_res = res;
351 
352  return global_res;
353 
354 }
355 
356 }
amrsolver::gauss_seidel
void gauss_seidel(int level, std::string name)
Definition: poisson.cpp:69
myOctree::scalar_fields
std::vector< std::string > scalar_fields
Definition: block.cpp:8
myOctree::Block::dy
double dy
Definition: block.h:57
myOctree::ny_block
int ny_block
Definition: block.cpp:11
myOctree::Field::name
std::string name
Definition: field.h:21
myOctree::Block::dx
double dx
Definition: block.h:57
myOctree::VecField
Template class for any vector field variable in the domain.
Definition: vecfield.h:13
myOctree::level_nodes
std::list< Octree * > level_nodes[20]
Definition: octreegrid.cpp:17
amrsolver::gauss_seidel_black
double gauss_seidel_black(Octree *node, Field *f, double global_res)
Definition: poisson.cpp:283
myOctree::Field
Template class for any scalar field variable in the domain.
Definition: field.h:11
amrsolver
Solver stuff.
Definition: amrsolver.cpp:15
myOctree::exchange_ghost_val
void exchange_ghost_val(int level, std::string name)
Definition: ghost.cpp:10
myOctree::pad
int pad
Definition: block.cpp:7
myOctree::Field::val
double *** val
Definition: field.h:20
myOctree::nz_block
int nz_block
Definition: block.cpp:12
myOctree::nx_block
int nx_block
Definition: block.cpp:10
amrsolver::gauss_seidel_red
double gauss_seidel_red(Octree *node, Field *f, double global_res)
Definition: poisson.cpp:207
myOctree::Octree
Class to store octree datastructure as nodes of the tree.
Definition: octree.h:26
myOctree::Octree::get_block_data
Block * get_block_data()
Definition: octree.cpp:20
myOctree::Block::dz
double dz
Definition: block.h:57
myOctree::create_lists_of_level_nodes
void create_lists_of_level_nodes()
Definition: octreegrid.cpp:43
amrsolver::jacobi
void jacobi(int level, std::string name)
Definition: poisson.cpp:23
amrsolver::jacobi_for_field
double jacobi_for_field(Octree *node, Field *f, double global_res)
Definition: poisson.cpp:130