Dendro  5.01
Dendro in Greek language means tree. The Dendro library is a large scale (262K cores on ORNL's Titan) distributed memory adaptive octree framework. The main goal of Dendro is to perform large scale multiphysics simulations efficeiently in mordern supercomputers. Dendro consists of efficient parallel data structures and algorithms to perform variational ( finite element) methods and finite difference mthods on 2:1 balanced arbitary adaptive octrees which enables the users to perform simulations raning from black holes (binary black hole mergers) to blood flow in human body, where applications ranging from relativity, astrophysics to biomedical engineering.
rkTransportUtils.h
1 //
2 // Created by milinda on 7/26/17.
8 //
9 
10 #ifndef SFCSORTBENCH_RKTRANSPORTUTILS_H
11 #define SFCSORTBENCH_RKTRANSPORTUTILS_H
12 
13 #include "mesh.h"
14 #include "block.h"
15 #include <iostream>
16 #include <vector>
17 #include "fdCoefficient.h"
18 
19 namespace adv_param
20 {
21  static const Point domain_min(-M_PI,-M_PI,-M_PI);
22  static const Point domain_max(M_PI,M_PI,M_PI);
23 }
24 
25 
33 template <typename T>
34 void grad(const ot::Mesh* pMesh,unsigned int dir,const T* uZipIn,T* uZipOut);
35 
36 
37 template <typename T>
38 void grad(const ot::Mesh* pMesh,unsigned int dir,const T* uZipIn,T* uZipOut)
39 {
40 
41 
42  const std::vector<ot::Block> blkList=pMesh->getLocalBlockList();
43  ot::TreeNode blkNode;
44  double h;
45  unsigned int blkNpe_1D;
46 
47  unsigned int grid_min=0;
48  unsigned int grid_max=(1u<<m_uiMaxDepth);
49  unsigned int paddWidth;
50  const unsigned int stencilSz=5;
51  unsigned int lx,ly,lz,offset;
52 
53 
54  if(dir==0)
55  { // compute the derivative in w.r.t x direction
56 
57 
58 
59  for(unsigned int blk=0;blk<blkList.size();blk++)
60  {
61 
62  blkNode=blkList[blk].getBlockNode();
63  paddWidth=blkList[blk].get1DPadWidth();
64  blkNpe_1D=blkList[blk].get1DArraySize();
65  h=1.0/(blkList[blk].computeDx(adv_param::domain_min,adv_param::domain_max));
66 
67  lx=blkList[blk].getAllocationSzX();
68  ly=blkList[blk].getAllocationSzY();
69  lz=blkList[blk].getAllocationSzZ();
70 
71  offset=blkList[blk].getOffset();
72 
73 
74  assert(blkNpe_1D>paddWidth);
75 
76  if(blkNode.minX()==grid_min)
77  { //std::cout<<"rank: "<<m_uiActiveRank<<" applying forward difference difference: "<<std::endl;
78 
79  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
80  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
81  {
82  uZipOut[offset+k*(ly*lx)+j*(lx)+paddWidth]=0;
83  for(unsigned int index=0;index<stencilSz;index++)
84  uZipOut[offset+k*(ly*lx)+j*(lx)+paddWidth]+=fd::D1_ORDER_4_FORWARD[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+paddWidth+index]*h;
85  }
86 
87 
88  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
89  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
90  for(unsigned int i=(paddWidth+1);i<(2*paddWidth);i++)
91  {
92  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
93  for(unsigned int index=0;index<stencilSz;index++)
94  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_UPWIND[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+i+index-1]*h;
95  }
96 
97 
98 
99  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
100  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
101  for(unsigned int i=2*paddWidth;i<(blkNpe_1D-2*paddWidth);i++)
102  {
103  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
104  for(unsigned int index=0;index<stencilSz;index++)
105  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+i+index-2]*h;
106  }
107 
108 
109 
110  if(blkNode.maxX()==grid_max)
111  {
112  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
113  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
114  for(unsigned int i=(blkNpe_1D-2*paddWidth);i<(blkNpe_1D-paddWidth-1);i++)
115  {
116  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
117  for(unsigned int index=0;index<stencilSz;index++)
118  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_DOWNWIND[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+i+index-3]*h;
119  }
120 
121 
122 
123  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
124  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
125  {
126  uZipOut[offset+k*(ly*lx)+j*(lx)+(blkNpe_1D-paddWidth-1)]=0;
127  for(unsigned int index=0;index<stencilSz;index++)
128  uZipOut[offset+k*(ly*lx)+j*(lx)+(blkNpe_1D-paddWidth-1)]+=fd::D1_ORDER_4_BACKWARD[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+(blkNpe_1D-paddWidth-1)+index-4]*h;
129  }
130 
131 
132 
133  }else
134  {
135  assert(blkNode.maxX()<grid_max);
136  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
137  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
138  for(unsigned int i=(blkNpe_1D-2*paddWidth);i<(blkNpe_1D-paddWidth);i++)
139  {
140  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
141  for(unsigned int index=0;index<stencilSz;index++)
142  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+i+index-2]*h;
143  }
144 
145 
146  }
147 
148  }else if(blkNode.maxX()==grid_max)
149  {
150 
151  assert(blkNode.minX()>grid_min);
152  assert((blkNpe_1D-2*paddWidth));
153  //std::cout<<"rank: "<<m_uiActiveRank<<" applying backward difference difference: "<<std::endl;
154  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
155  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
156  for(unsigned int i=paddWidth;i<(blkNpe_1D-2*paddWidth);i++)
157  {
158  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
159  for(unsigned int index=0;index<stencilSz;index++)
160  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+i+index-2]*h;
161 
162  }
163 
164 
165  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
166  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
167  for(unsigned int i=(blkNpe_1D-2*paddWidth);i<(blkNpe_1D-paddWidth-1);i++)
168  {
169  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
170  for(unsigned int index=0;index<stencilSz;index++)
171  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_DOWNWIND[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+i+index-3]*h;
172  }
173 
174 
175 
176  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
177  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
178  {
179  uZipOut[offset+k*(ly*lx)+j*(lx)+(blkNpe_1D-paddWidth-1)]=0;
180  for(unsigned int index=0;index<stencilSz;index++)
181  uZipOut[offset+k*(ly*lx)+j*(lx)+(blkNpe_1D-paddWidth-1)]+=fd::D1_ORDER_4_BACKWARD[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+(blkNpe_1D-paddWidth-1)+index-4]*h;
182  }
183 
184 
185 
186  }else
187  {
188  assert(blkNode.minX()>grid_min && blkNode.maxX()<grid_max);
189 
190  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
191  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
192  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
193  {
194  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
195  for(unsigned int index=0;index<stencilSz;index++)
196  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+k*(ly*lx)+j*(lx)+(i+index-2)]*h;
197  }
198 
199 
200  }
201 
202 
203  }
204  }else if(dir==1)
205  { // compute the derivative in w.r.t y direction
206 
207  for(unsigned int blk=0;blk<blkList.size();blk++)
208  {
209 
210  blkNode=blkList[blk].getBlockNode();
211  paddWidth=blkList[blk].get1DPadWidth();
212  blkNpe_1D=blkList[blk].get1DArraySize();
213  h=1.0/(blkList[blk].computeDy(adv_param::domain_min,adv_param::domain_max));
214 
215  lx=blkList[blk].getAllocationSzX();
216  ly=blkList[blk].getAllocationSzY();
217  lz=blkList[blk].getAllocationSzZ();
218 
219  offset=blkList[blk].getOffset();
220 
221 
222  assert(blkNpe_1D>paddWidth);
223 
224  if(blkNode.minY()==grid_min)
225  { //std::cout<<"rank: "<<m_uiActiveRank<<" applying forward difference difference: "<<std::endl;
226 
227  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
228  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
229  {
230  uZipOut[offset+k*(ly*lx)+(paddWidth)*(lx)+i]=0;
231  for(unsigned int index=0;index<stencilSz;index++)
232  uZipOut[offset+k*(ly*lx)+(paddWidth)*(lx)+i]+=fd::D1_ORDER_4_FORWARD[index]*uZipIn[offset+k*(ly*lx)+(paddWidth+index)*(lx)+i]*h;
233  }
234 
235 
236  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
237  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
238  for(unsigned int j=(paddWidth+1);j<(2*paddWidth);j++)
239  {
240  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
241  for(unsigned int index=0;index<stencilSz;index++)
242  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_UPWIND[index]*uZipIn[offset+k*(ly*lx)+(j+index-1)*(lx)+i]*h;
243  }
244 
245 
246 
247  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
248  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
249  for(unsigned int j=2*paddWidth;j<(blkNpe_1D-2*paddWidth);j++)
250  {
251  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
252  for(unsigned int index=0;index<stencilSz;index++)
253  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+k*(ly*lx)+(j+index-2)*(lx)+i]*h;
254  }
255 
256 
257 
258  if(blkNode.maxY()==grid_max)
259  {
260  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
261  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
262  for(unsigned int j=(blkNpe_1D-2*paddWidth);j<(blkNpe_1D-paddWidth-1);j++)
263  {
264  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
265  for(unsigned int index=0;index<stencilSz;index++)
266  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_DOWNWIND[index]*uZipIn[offset+k*(ly*lx)+(j+index-3)*(lx)+i]*h;
267  }
268 
269 
270 
271  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
272  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
273  {
274  uZipOut[offset+k*(ly*lx)+(blkNpe_1D-paddWidth-1)*(lx)+i]=0;
275  for(unsigned int index=0;index<stencilSz;index++)
276  uZipOut[offset+k*(ly*lx)+(blkNpe_1D-paddWidth-1)*(lx)+i]+=fd::D1_ORDER_4_BACKWARD[index]*uZipIn[offset+k*(ly*lx)+((blkNpe_1D-paddWidth-1)+index-4)*(lx)+i]*h;
277  }
278 
279 
280 
281  }else
282  {
283  assert(blkNode.maxY()<grid_max);
284  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
285  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
286  for(unsigned int j=(blkNpe_1D-2*paddWidth);j<(blkNpe_1D-paddWidth);j++)
287  {
288  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
289  for(unsigned int index=0;index<stencilSz;index++)
290  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+k*(ly*lx)+(j+index-2)*(lx)+i]*h;
291  }
292 
293 
294  }
295 
296  }else if(blkNode.maxY()==grid_max)
297  {
298 
299  assert(blkNode.minY()>grid_min);
300  assert((blkNpe_1D-2*paddWidth));
301  //std::cout<<"rank: "<<m_uiActiveRank<<" applying backward difference difference: "<<std::endl;
302  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
303  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
304  for(unsigned int j=paddWidth;j<(blkNpe_1D-2*paddWidth);j++)
305  {
306  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
307  for(unsigned int index=0;index<stencilSz;index++)
308  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+k*(ly*lx)+(j+index-2)*(lx)+i]*h;
309 
310  }
311 
312 
313  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
314  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
315  for(unsigned int j=(blkNpe_1D-2*paddWidth);j<(blkNpe_1D-paddWidth-1);j++)
316  {
317  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
318  for(unsigned int index=0;index<stencilSz;index++)
319  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_DOWNWIND[index]*uZipIn[offset+k*(ly*lx)+(j+index-3)*(lx)+i]*h;
320  }
321 
322 
323 
324  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
325  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
326  {
327  uZipOut[offset+k*(ly*lx)+(blkNpe_1D-paddWidth-1)*(lx)+i]=0;
328  for(unsigned int index=0;index<stencilSz;index++)
329  uZipOut[offset+k*(ly*lx)+(blkNpe_1D-paddWidth-1)*(lx)+i]+=fd::D1_ORDER_4_BACKWARD[index]*uZipIn[offset+k*(ly*lx)+((blkNpe_1D-paddWidth-1)+index-4)*(lx)+i]*h;
330  }
331 
332 
333 
334  }else
335  {
336  assert(blkNode.minY()>grid_min && blkNode.maxY()<grid_max);
337 
338  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
339  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
340  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
341  {
342  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
343  for(unsigned int index=0;index<stencilSz;index++)
344  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+k*(ly*lx)+(j+index-2)*(lx)+i]*h;
345  }
346 
347 
348  }
349 
350 
351  }
352 
353  }else if(dir==2)
354  { // compute the derivative in w.r.t z direction
355 
356  for(unsigned int blk=0;blk<blkList.size();blk++)
357  {
358 
359  blkNode=blkList[blk].getBlockNode();
360  paddWidth=blkList[blk].get1DPadWidth();
361  blkNpe_1D=blkList[blk].get1DArraySize();
362  h=1.0/(blkList[blk].computeDz(adv_param::domain_min,adv_param::domain_max));
363 
364  lx=blkList[blk].getAllocationSzX();
365  ly=blkList[blk].getAllocationSzY();
366  lz=blkList[blk].getAllocationSzZ();
367 
368  offset=blkList[blk].getOffset();
369 
370 
371  assert(blkNpe_1D>paddWidth);
372 
373  if(blkNode.minZ()==grid_min)
374  { //std::cout<<"rank: "<<m_uiActiveRank<<" applying forward difference difference: "<<std::endl;
375 
376  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
377  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
378  {
379  uZipOut[offset+(paddWidth)*(ly*lx)+(j)*(lx)+i]=0;
380  for(unsigned int index=0;index<stencilSz;index++)
381  uZipOut[offset+(paddWidth)*(ly*lx)+(j)*(lx)+i]+=fd::D1_ORDER_4_FORWARD[index]*uZipIn[offset+(paddWidth+index)*(ly*lx)+j*(lx)+i]*h;
382  }
383 
384 
385  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
386  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
387  for(unsigned int k=(paddWidth+1);k<(2*paddWidth);k++)
388  {
389  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
390  for(unsigned int index=0;index<stencilSz;index++)
391  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_UPWIND[index]*uZipIn[offset+(k+index-1)*(ly*lx)+j*(lx)+i]*h;
392  }
393 
394 
395 
396  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
397  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
398  for(unsigned int k=2*paddWidth;k<(blkNpe_1D-2*paddWidth);k++)
399  {
400  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
401  for(unsigned int index=0;index<stencilSz;index++)
402  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+(k+index-2)*(ly*lx)+j*(lx)+i]*h;
403  }
404 
405 
406 
407  if(blkNode.maxZ()==grid_max)
408  {
409  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
410  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
411  for(unsigned int k=(blkNpe_1D-2*paddWidth);k<(blkNpe_1D-paddWidth-1);k++)
412  {
413  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
414  for(unsigned int index=0;index<stencilSz;index++)
415  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_DOWNWIND[index]*uZipIn[offset+(k+index-3)*(ly*lx)+j*(lx)+i]*h;
416  }
417 
418 
419 
420  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
421  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
422  {
423  uZipOut[offset+(blkNpe_1D-paddWidth-1)*(ly*lx)+j*(lx)+i]=0;
424  for(unsigned int index=0;index<stencilSz;index++)
425  uZipOut[offset+(blkNpe_1D-paddWidth-1)*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_BACKWARD[index]*uZipIn[offset+((blkNpe_1D-paddWidth-1)+index-4)*(ly*lx)+j*(lx)+i]*h;
426  }
427 
428 
429 
430  }else
431  {
432  assert(blkNode.maxY()<grid_max);
433  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
434  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
435  for(unsigned int k=(blkNpe_1D-2*paddWidth);k<(blkNpe_1D-paddWidth);k++)
436  {
437  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
438  for(unsigned int index=0;index<stencilSz;index++)
439  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+(k+index-2)*(ly*lx)+j*(lx)+i]*h;
440  }
441 
442 
443  }
444 
445  }else if(blkNode.maxZ()==grid_max)
446  {
447 
448  assert(blkNode.minZ()>grid_min);
449  assert((blkNpe_1D-2*paddWidth));
450  //std::cout<<"rank: "<<m_uiActiveRank<<" applying backward difference difference: "<<std::endl;
451  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
452  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
453  for(unsigned int k=paddWidth;k<(blkNpe_1D-2*paddWidth);k++)
454  {
455  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
456  for(unsigned int index=0;index<stencilSz;index++)
457  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+(k+index-2)*(ly*lx)+j*(lx)+i]*h;
458 
459  }
460 
461 
462  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
463  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
464  for(unsigned int k=(blkNpe_1D-2*paddWidth);k<(blkNpe_1D-paddWidth-1);k++)
465  {
466  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
467  for(unsigned int index=0;index<stencilSz;index++)
468  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_DOWNWIND[index]*uZipIn[offset+(k+index-3)*(ly*lx)+j*(lx)+i]*h;
469  }
470 
471 
472 
473  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
474  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
475  {
476  uZipOut[offset+(blkNpe_1D-paddWidth-1)*(ly*lx)+j*(lx)+i]=0;
477  for(unsigned int index=0;index<stencilSz;index++)
478  uZipOut[offset+(blkNpe_1D-paddWidth-1)*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_BACKWARD[index]*uZipIn[offset+((blkNpe_1D-paddWidth-1)+index-4)*(ly*lx)+j*(lx)+i]*h;
479  }
480 
481 
482 
483  }else
484  {
485  assert(blkNode.minZ()>grid_min && blkNode.maxZ()<grid_max);
486 
487  for(unsigned int j=paddWidth;j<(blkNpe_1D-paddWidth);j++)
488  for(unsigned int i=paddWidth;i<(blkNpe_1D-paddWidth);i++)
489  for(unsigned int k=paddWidth;k<(blkNpe_1D-paddWidth);k++)
490  {
491  uZipOut[offset+k*(ly*lx)+j*(lx)+i]=0;
492  for(unsigned int index=0;index<stencilSz;index++)
493  uZipOut[offset+k*(ly*lx)+j*(lx)+i]+=fd::D1_ORDER_4_CENTERED[index]*uZipIn[offset+(k+index-2)*(ly*lx)+j*(lx)+i]*h;
494  }
495 
496 
497  }
498 
499 
500  }
501 
502 
503 
504  }else
505  {
506  std::cout<<" unknown stencil direction "<<std::endl;
507  }
508 
509 
510 
511 }
512 
513 
514 
515 #endif //SFCSORTBENCH_RKTRANSPORTUTILS_H
ot::TreeNode::minX
unsigned int minX() const
returns true min corner(anchor) and the max corner coordinates of a octant.
Definition: TreeNode.h:164
ot::TreeNode
A class to manage octants.
Definition: TreeNode.h:35
ot::Mesh
Definition: mesh.h:179
Point
A point class.
Definition: point.h:36
adv_param
Contains gradient computation functions for the advection equation.
Definition: rkTransportUtils.h:19
ot::Mesh::getLocalBlockList
const std::vector< ot::Block > & getLocalBlockList() const
returns const list of local blocks (regular grids) for the consdering mesh.
Definition: mesh.h:1084
Generated by   doxygen 1.8.13