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.
operators.h
1 //
2 // Created by milinda on 12/26/17.
12 //
13 
14 #ifndef SFCSORTBENCH_OPERATORS_H
15 #define SFCSORTBENCH_OPERATORS_H
16 
17 #include "TreeNode.h"
18 #include "refel.h"
19 #include <functional>
20 #include "tensor.h"
21 #include "workspace.h"
22 #include "mathUtils.h"
23 #include "dendroProfileParams.h"
24 
25 #define GRIDX_TO_X(xg) (((Rx/RgX)*(xg-bssn::BSSN_OCTREE_MIN[0]))+bssn::BSSN_COMPD_MIN[0])
26 #define GRIDY_TO_Y(yg) (((Ry/RgY)*(yg-bssn::BSSN_OCTREE_MIN[1]))+bssn::BSSN_COMPD_MIN[1])
27 #define GRIDZ_TO_Z(zg) (((Rz/RgZ)*(zg-bssn::BSSN_OCTREE_MIN[2]))+bssn::BSSN_COMPD_MIN[2])
28 
29 #define X_TO_GRIDX(xc) (((RgX/Rx)*(xc-bssn::BSSN_COMPD_MIN[0]))+bssn::BSSN_OCTREE_MIN[0])
30 #define Y_TO_GRIDY(yc) (((RgY/Ry)*(yc-bssn::BSSN_COMPD_MIN[1]))+bssn::BSSN_OCTREE_MIN[1])
31 #define Z_TO_GRIDZ(zc) (((RgZ/Rz)*(zc-bssn::BSSN_COMPD_MIN[2]))+bssn::BSSN_OCTREE_MIN[2])
32 
33 
34 
35 
36 namespace fem
37 {
38  namespace operators
39  {
40 
41  namespace poisson
42  {
43 
54  template <typename T>
55  void elementMassMatvec(const T* f_rhs, const ot::TreeNode & elem, const RefElement * refEl, T* out)
56  {
57 
58 #ifdef MATVEC_PROFILE
59  dendro::timer::sfcmatvec::t_elemMvec.start();
60 #endif
61 
62  const double * Q1d=refEl->getQ1d();
63  const double * QT1d=refEl->getQT1d();
64  const double * Dg=refEl->getDg1d();
65  const double * W1d=refEl->getWgq();
66 
67  const unsigned int eleOrder=refEl->getOrder();
68  const unsigned int nPe=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
69  const unsigned int nrp=eleOrder+1;
70 
71  const double szX=fem::domain::gridX_to_X(elem.maxX())-fem::domain::gridX_to_X(elem.minX());
72  const double szY=fem::domain::gridY_to_Y(elem.maxY())-fem::domain::gridY_to_Y(elem.minY());
73  const double szZ=fem::domain::gridZ_to_Z(elem.maxZ())-fem::domain::gridZ_to_Z(elem.minZ());
74 
75 
76  const double refElSz=refEl->getElementSz();
77 
78  // interpolate to quadrature points.
79  DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,Q1d,f_rhs,imV1);
80  DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,Q1d,imV1,imV2);
81  DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,Q1d,imV2,out);
82 
83 
84 
85  const double Jx = 1.0/(refElSz/(double (szX)));
86  const double Jy = 1.0/(refElSz/(double (szY)));
87  const double Jz = 1.0/(refElSz/(double (szZ)));
88 
89  //std::cout<<"Mass: elem: "<<elem<<" ele Sz: "<<(elem.maxX()-elem.minX())<<" szX: "<<szX<<" Jx: "<<Jx<<" J: "<<(Jx*Jy*Jz)<<std::endl;
90 
91  for(unsigned int k=0;k<(eleOrder+1);k++)
92  for(unsigned int j=0;j<(eleOrder+1);j++)
93  for(unsigned int i=0;i<(eleOrder+1);i++)
94  out[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]*=(Jx*Jy*Jz*W1d[i]*W1d[j]*W1d[k]);
95 
96 
97  // apply transpose operator
98  DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,QT1d,out,imV1);
99  DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,QT1d,imV1,imV2);
100  DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,QT1d,imV2,out);
101 
102 #ifdef FEM_ACCUMILATE_ONES_TEST
103  for(unsigned int k=0;k<(eleOrder+1);k++)
104  for(unsigned int j=0;j<(eleOrder+1);j++)
105  for(unsigned int i=0;i<(eleOrder+1);i++)
106  out[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]=1.0;
107 
108 #endif
109  /*if(elem.minX()==0 && elem.minY()==0 && elem.minZ()==0)
110  {
111  std::cout<<"elem: "<<elem<<std::endl;
112  for(unsigned int k=0;k<(eleOrder+1);k++)
113  for(unsigned int j=0;j<(eleOrder+1);j++)
114  for(unsigned int i=0;i<(eleOrder+1);i++)
115  std::cout<<"i: "<<i<<" j: "<<j<<" k: "<<k<<" MF "<<out[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]<<" F "<<f_rhs[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]<<std::endl;
116 
117  }*/
118 
119 #ifdef MATVEC_PROFILE
120  dendro::timer::sfcmatvec::t_elemMvec.stop();
121 #endif
122 
123  }
124 
129  template <typename T>
130  void elementStiffnessMatvec(const T* in, const ot::TreeNode & elem, const RefElement * refEl, T* out)
131  {
132 
133 #ifdef MATVEC_PROFILE
134  dendro::timer::sfcmatvec::t_elemMvec.start();
135 #endif
136 
137  const double * Q1d=refEl->getQ1d();
138  const double * QT1d=refEl->getQT1d();
139  const double * Dg=refEl->getDg1d();
140  const double * DgT=refEl->getDgT1d();
141  const double * W1d=refEl->getWgq();
142 
143  const unsigned int eleOrder=refEl->getOrder();
144  const unsigned int nPe=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
145  const unsigned int nrp=eleOrder+1;
146 
147  const unsigned int sz=1u<<(m_uiMaxDepth-elem.getLevel()); // curent element size
148  const double refElSz=refEl->getElementSz();
149  //x derivative
150  DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,Dg,in,imV1);
151  DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,Q1d,imV1,imV2);
152  DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,Q1d,imV2,Qx);
153 
154  //y derivative
155  DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,Q1d,in,imV1);
156  DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,Dg,imV1,imV2);
157  DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,Q1d,imV2,Qy);
158 
159  //z derivative
160  DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,Q1d,in,imV1);
161  DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,Q1d,imV1,imV2);
162  DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,Dg,imV2,Qz);
163 
164 
165  const double szX=fem::domain::gridX_to_X(elem.maxX())-fem::domain::gridX_to_X(elem.minX());
166  const double szY=fem::domain::gridY_to_Y(elem.maxY())-fem::domain::gridY_to_Y(elem.minY());
167  const double szZ=fem::domain::gridZ_to_Z(elem.maxZ())-fem::domain::gridZ_to_Z(elem.minZ());
168 
169  const double Jx = 1.0/(refElSz/(double (szX)));
170  const double Jy = 1.0/(refElSz/(double (szY)));
171  const double Jz = 1.0/(refElSz/(double (szZ)));
172 
173  //std::cout<<"Stifness: elem: "<<elem<<" ele Sz: "<<(elem.maxX()-elem.minX())<<" szX: "<<szX<<" Jx: "<<Jx<<" J: "<<(Jx*Jy*Jz)<<std::endl;
174 
175  for(unsigned int k=0;k<(eleOrder+1);k++)
176  for(unsigned int j=0;j<(eleOrder+1);j++)
177  for(unsigned int i=0;i<(eleOrder+1);i++)
178  {
179  Qx[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]*=( ((Jy*Jz)/Jx)*W1d[i]*W1d[j]*W1d[k]);
180  Qy[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]*=( ((Jx*Jz)/Jy)*W1d[i]*W1d[j]*W1d[k]);
181  Qz[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]*=( ((Jx*Jy)/Jz)*W1d[i]*W1d[j]*W1d[k]);
182  }
183 
184 
185 
186  DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,DgT,Qx,imV1);
187  DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,QT1d,imV1,imV2);
188  DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,QT1d,imV2,Qx);
189 
190  DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,QT1d,Qy,imV1);
191  DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,DgT,imV1,imV2);
192  DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,QT1d,imV2,Qy);
193 
194  DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,QT1d,Qz,imV1);
195  DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,QT1d,imV1,imV2);
196  DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,DgT,imV2,Qz);
197 
198  for(unsigned int i=0;i<nPe;i++)
199  out[i]=Qx[i]+Qy[i]+Qz[i];
200 
201 
202 #ifdef MATVEC_PROFILE
203  dendro::timer::sfcmatvec::t_elemMvec.stop();
204 #endif
205 
206 
207  }
208 
209 
210 
211  template <typename T>
212  void applyDirichletBoundary(const ot::Mesh* mesh,T* in,T bdyValue)
213  {
214 
215  const unsigned int nPe=mesh->getNumNodesPerElement();
216  const unsigned int eleOrder=mesh->getElementOrder();
217  const ot::TreeNode* pNodes = &(*(mesh->getAllElements().begin()));
218 
219  const unsigned int * e2n=&(*(mesh->getE2NMapping().begin()));
220  const unsigned int * e2n_dg=&(*(mesh->getE2NMapping_DG().begin()));
221 
222  const unsigned int d_min=0;
223  const unsigned int d_max=(1u<<m_uiMaxDepth);
224 
225 
226  // apply zero-dirichlet bdy conditions here.
227  const unsigned int nx=(eleOrder+1);
228  const unsigned int ny=(eleOrder+1)*(eleOrder+1);
229  const unsigned int nz=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
230  const double bdy_value=bdyValue;
231  double x,y,z;
232 
233  for(unsigned int ele=mesh->getElementLocalBegin();ele<mesh->getElementLocalEnd();ele++)
234  {
235  x=pNodes[ele].minX();
236  y=pNodes[ele].minY();
237  z=pNodes[ele].minZ();
238 
239  if(x==d_min)
240  {
241  for(unsigned int k=0;k<(eleOrder+1);k++)
242  for(unsigned int j=0;j<(eleOrder+1);j++)
243  {
244  if( (!mesh->isNodeHanging(ele,0,j,k)) && (mesh->isNodeLocal(ele,0,j,k)))
245  in[e2n[ele*nPe+k*ny+j*nx+0]]=bdy_value;
246  }
247  }
248 
249  if(y==d_min)
250  {
251  for(unsigned int k=0;k<(eleOrder+1);k++)
252  for(unsigned int i=0;i<(eleOrder+1);i++)
253  {
254  if( (!mesh->isNodeHanging(ele,i,0,k)) && (mesh->isNodeLocal(ele,i,0,k)))
255  in[e2n[ele*nPe+k*ny+0*nx+i]]=bdy_value;
256  }
257  }
258 
259 
260  if(z==d_min)
261  {
262  for(unsigned int j=0;j<(eleOrder+1);j++)
263  for(unsigned int i=0;i<(eleOrder+1);i++)
264  {
265  if( (!mesh->isNodeHanging(ele,i,j,0)) && (mesh->isNodeLocal(ele,i,j,0)))
266  in[e2n[ele*nPe+0*ny+j*nx+i]]=bdy_value;
267  }
268  }
269 
270  x=pNodes[ele].maxX();
271  y=pNodes[ele].maxY();
272  z=pNodes[ele].maxZ();
273 
274  if(x==d_max)
275  {
276  for(unsigned int k=0;k<(eleOrder+1);k++)
277  for(unsigned int j=0;j<(eleOrder+1);j++)
278  {
279  if( (!mesh->isNodeHanging(ele,eleOrder,j,k)) && (mesh->isNodeLocal(ele,eleOrder,j,k)))
280  in[e2n[ele*nPe+k*ny+j*nx+eleOrder]]=bdy_value;
281  }
282  }
283 
284  if(y==d_max)
285  {
286  for(unsigned int k=0;k<(eleOrder+1);k++)
287  for(unsigned int i=0;i<(eleOrder+1);i++)
288  {
289  if( (!mesh->isNodeHanging(ele,i,eleOrder,k)) && (mesh->isNodeLocal(ele,i,eleOrder,k)))
290  in[e2n[ele*nPe+k*ny+eleOrder*nx+i]]=bdy_value;
291  }
292  }
293 
294 
295  if(z==d_max)
296  {
297  for(unsigned int j=0;j<(eleOrder+1);j++)
298  for(unsigned int i=0;i<(eleOrder+1);i++)
299  {
300  if( (!mesh->isNodeHanging(ele,i,j,eleOrder)) && (mesh->isNodeLocal(ele,i,j,eleOrder)))
301  in[e2n[ele*nPe+eleOrder*ny+j*nx+i]]=bdy_value;
302  }
303 
304  }
305 
306 
307  }
308  }
309 
310 
311 
312 
320  template <typename T>
321  void matvec(ot::Mesh* mesh,T* in, T* out,bool applyBoundary=true)
322  {
323  if(mesh->isActive())
324  {
325 
326  const std::vector<unsigned int> level1Ghost=mesh->getLevel1GhostElementIndices();
327 
328  mesh->performGhostExchange(in);
329 
330  const unsigned int nPe=mesh->getNumNodesPerElement();
331  const unsigned int eleOrder=mesh->getElementOrder();
332 
333 
334  T* elementalVec=new T[nPe];
335  T* elementalMVec=new T[nPe];
336 
337  const unsigned int nx=(eleOrder+1);
338  const unsigned int ny=(eleOrder+1)*(eleOrder+1);
339  const unsigned int nz=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
340  const ot::TreeNode* pNodes = &(*(mesh->getAllElements().begin()));
341  const RefElement *refEl=mesh->getReferenceElement();
342 
343  unsigned int nodeLookup,ownerID,ii_x,jj_y,kk_z;
344  for(unsigned int node=mesh->getNodeLocalBegin();node<mesh->getNodeLocalEnd();node++)
345  out[node]=(T)0;
346 
347  mesh->performGhostExchange(out);
348 
349  // iterate over ghost elements . (later we can use this to potentially overlap computation & communication)
350  for(unsigned int gIndex=0;gIndex<level1Ghost.size();gIndex++)
351  {
352  unsigned int ele=level1Ghost[gIndex];
353  mesh->getElementNodalValues(in,elementalVec,ele);
354  operators::poisson::elementStiffnessMatvec(elementalVec,pNodes[ele],refEl,elementalMVec);
355  mesh->computeElementalContribution(elementalMVec,out,ele);
356  }
357 
358 
359 
360  // iterate over local elements.
361  for(unsigned int ele=mesh->getElementLocalBegin();ele<mesh->getElementLocalEnd();ele++)
362  {
363  mesh->getElementNodalValues(in,elementalVec,ele);
364 
365  /* if(ele==0)
366  for(unsigned int node=0;node<125;node++)
367  std::cout<<" node value in: "<<elementalVec[node]<<std::endl;*/
368 
369  operators::poisson::elementStiffnessMatvec(elementalVec,pNodes[ele],refEl,elementalMVec);
370 
371  /* if(ele==0)
372  for(unsigned int node=0;node<125;node++)
373  std::cout<<" node value out: "<<elementalMVec[node]<<std::endl;*/
374 
375  mesh->computeElementalContribution(elementalMVec,out,ele);
376 
377  }
378 
379  delete [] elementalVec;
380  delete [] elementalMVec;
381  mesh->performGhostExchange(out);
382 
383 
384  }
385  }
386 
387 
395  template <typename T>
396  void computeRHS(ot::Mesh* mesh,T* in, T* out)
397  {
398 
399 
400  if(mesh->isActive())
401  {
402  //applyDirichletBoundary(mesh,in,0.0);
403  mesh->performGhostExchange(in);
404  const unsigned int nPe=mesh->getNumNodesPerElement();
405  const unsigned int eleOrder=mesh->getElementOrder();
406  const std::vector<unsigned int> level1Ghost=mesh->getLevel1GhostElementIndices();
407 
408 
409  T* elementalVec=new T[nPe];
410  T* elementalMVec=new T[nPe];
411 
412  const unsigned int nx=(eleOrder+1);
413  const unsigned int ny=(eleOrder+1)*(eleOrder+1);
414  const unsigned int nz=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
415  const ot::TreeNode* pNodes = &(*(mesh->getAllElements().begin()));
416  const RefElement *refEl=mesh->getReferenceElement();
417 
418  /*const double * W1d=refEl->getWgq();
419  for(unsigned int i=0;i<5;i++)
420  std::cout<<" w["<<i<<"]: "<<W1d[i]<<std::endl;*/
421 
422  const unsigned int * e2n=&(*(mesh->getE2NMapping().begin()));
423  const unsigned int * e2n_dg=&(*(mesh->getE2NMapping_DG().begin()));
424  unsigned int nodeLookup,ownerID,ii_x,jj_y,kk_z;
425 
426 
427  for(unsigned int node=mesh->getNodeLocalBegin();node<mesh->getNodeLocalEnd();node++)
428  out[node]=(T)0.0;
429 
430  mesh->performGhostExchange(out);
431 
432  // iterate over ghost elements . (later we can use this to potentially overlap computation & communication)
433  for(unsigned int gIndex=0;gIndex<level1Ghost.size();gIndex++)
434  {
435  unsigned int ele=level1Ghost[gIndex];
436 
437 #ifdef MATVEC_PROFILE
438  dendro::timer::sfcmatvec::t_p2cInterp.start();
439 #endif
440  mesh->getElementNodalValues(in,elementalVec,ele);
441 #ifdef MATVEC_PROFILE
442  dendro::timer::sfcmatvec::t_p2cInterp.stop();
443 #endif
444 
445  operators::poisson::elementMassMatvec(elementalVec,pNodes[ele],refEl,elementalMVec);
446 
447 
448 #ifdef MATVEC_PROFILE
449  dendro::timer::sfcmatvec::t_c2pInterp.start();
450 #endif
451 
452  mesh->computeElementalContribution(elementalMVec,out,ele);
453 
454 #ifdef MATVEC_PROFILE
455  dendro::timer::sfcmatvec::t_c2pInterp.stop();
456 #endif
457 
458  }
459 
460  // iterate over local elements
461  for(unsigned int ele=mesh->getElementLocalBegin();ele<mesh->getElementLocalEnd();ele++)
462  {
463 
464 #ifdef MATVEC_PROFILE
465  dendro::timer::sfcmatvec::t_p2cInterp.start();
466 #endif
467  mesh->getElementNodalValues(in,elementalVec,ele);
468 
469 #ifdef MATVEC_PROFILE
470  dendro::timer::sfcmatvec::t_p2cInterp.stop();
471 #endif
472  operators::poisson::elementMassMatvec(elementalVec,pNodes[ele],refEl,elementalMVec);
473 
474 
475 #ifdef MATVEC_PROFILE
476  dendro::timer::sfcmatvec::t_c2pInterp.start();
477 #endif
478  mesh->computeElementalContribution(elementalMVec,out,ele);
479 
480 #ifdef MATVEC_PROFILE
481  dendro::timer::sfcmatvec::t_c2pInterp.stop();
482 #endif
483 
484 
485  }
486 
487  mesh->performGhostExchange(out);
488 
489  delete [] elementalVec;
490  delete [] elementalMVec;
491 
492 
493  }
494 
495  }
496 
497 
498 
499  } // end of namespace possoin
500 
501  } // end of namespace operators
502 
503 
504 } // end of namespace fem
505 
506 
507 
508 #endif //SFCSORTBENCH_OPERATORS_H
const std::vector< unsigned int > & getLevel1GhostElementIndices() const
returns the Level 1 ghost element indices.
Definition: mesh.h:1061
void computeElementalContribution(const T *in, T *out, unsigned int elementID) const
: Computes the contribution of elemental nodal values to the parent elements if it is hanging...
unsigned int minX() const
returns true min corner(anchor) and the max corner coordinates of a octant.
Definition: TreeNode.h:164
A class to manage octants.
Definition: TreeNode.h:35
void performGhostExchange(std::vector< T > &vec)
Perform the ghost exchange for the vector vec.
Definition: mesh.h:179
const std::vector< unsigned int > & getE2NMapping_DG() const
Definition: mesh.h:1082
unsigned int getNodeLocalBegin() const
return the location of local node begin
Definition: mesh.h:1043
bool isNodeLocal(unsigned int eleID, unsigned int ix, unsigned int jy, unsigned int kz) const
: Returns true if the specified node (e,i,j,k) is local.
Definition: mesh.h:1358
workspace variables allocations
Definition: matvec.h:44
unsigned int getElementLocalBegin() const
return the begin location of element local
Definition: mesh.h:1030
void getElementNodalValues(const T *vec, T *nodalValues, unsigned int elementID) const
: Returns the nodal values of a given element for a given variable vector.
const std::vector< ot::TreeNode > & getAllElements() const
returns the pointer to All elements array.
Definition: mesh.h:1058
unsigned int getNumNodesPerElement() const
return the number of nodes per element.
Definition: mesh.h:1090
const RefElement * getReferenceElement() const
returns const pointer to reference element
Definition: mesh.h:1114
unsigned int getLevel() const
return the level of the of the octant
Definition: TreeNode.h:387
bool isNodeHanging(unsigned int eleID, unsigned int ix, unsigned int jy, unsigned int kz) const
: Returns true if the specified node (e,i,j,k) is hanging.
Definition: mesh.cpp:8943
const std::vector< unsigned int > & getE2NMapping() const
Definition: mesh.h:1073
bool isActive() const
returns true if mesh is active
Definition: mesh.h:1173
unsigned int getNodeLocalEnd() const
return the location of local node end
Definition: mesh.h:1045
unsigned int getElementOrder() const
returns the order of an element
Definition: mesh.h:1093
unsigned int getElementLocalEnd() const
return the end location of element local
Definition: mesh.h:1032
Definition: refel.h:69