diff options
author | Erik Schnetter <schnetter@cct.lsu.edu> | 2010-04-27 10:21:08 -0500 |
---|---|---|
committer | Erik Schnetter <schnetter@cct.lsu.edu> | 2010-04-27 10:21:08 -0500 |
commit | 010cc61a2d2a6e1b8cb5dda17ac6b11eaf1ac0da (patch) | |
tree | 695d19c112b4607d29a8dec9d38793cc7ad111ba /Carpet/CarpetReduce | |
parent | 5c9f5be27278e3b5e12ad5f795f3d66b78e7296a (diff) |
CarpetReduce: Set up mask from new dh class entries
Simplify mask setup code and improve its performance by using the new
dh classes and their new entries.
Diffstat (limited to 'Carpet/CarpetReduce')
-rw-r--r-- | Carpet/CarpetReduce/src/mask_carpet.cc | 257 |
1 files changed, 73 insertions, 184 deletions
diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc index 260afb23d..e469d105b 100644 --- a/Carpet/CarpetReduce/src/mask_carpet.cc +++ b/Carpet/CarpetReduce/src/mask_carpet.cc @@ -22,28 +22,6 @@ namespace CarpetMask { - void ibbox2iminimax (ibbox const& ext, // component extent - ibbox const& box, // this bbox - ivect const& lsh, // component's lsh - ivect& imin, ivect& imax) - { - ivect const izero = ivect(0); - - assert (all ((box.lower() - ext.lower() ) >= 0)); - assert (all ((box.upper() - ext.lower() + ext.stride()) >= 0)); - assert (all ((box.lower() - ext.lower() ) % ext.stride() == 0)); - assert (all ((box.upper() - ext.lower() + ext.stride()) % ext.stride() == 0)); - - imin = (box.lower() - ext.lower() ) / ext.stride(); - imax = (box.upper() - ext.lower() + ext.stride()) / ext.stride(); - - assert (all (izero <= imin)); - assert (box.empty() xor all (imin <= imax)); - assert (all (imax <= lsh)); - } - - - /** * Reduce the weight on the current and the next coarser level to * make things consistent. Set the weight to 0 inside the @@ -76,86 +54,27 @@ namespace CarpetMask { - // Calculate the union of all refined regions - ibset active; - for (int c=0; c<hh.components(reflevel); ++c) { - // refined |= hh.extent(mglevel,reflevel,c); - ibset this_active; - dh::dboxes const& box = dd.boxes.AT(mglevel).AT(reflevel).AT(c); - dh::dboxes::ibboxs2ibset (box.active, box.numactive, this_active); - active += this_active; - } - active.normalize(); - - // Calculate the union of all coarse regions - ibset parent; - for (int c=0; c<hh.components(reflevel-1); ++c) { - parent |= hh.extent(mglevel,reflevel-1,c) - .expand(ivect(reffact-1),ivect(reffact-1)).contracted_for(base); - } - parent.normalize(); - - // Subtract the refined region - ibset notrefined = parent - active; - notrefined.normalize(); - - // Enlarge this set - ibset enlarged[dim]; - for (int d=0; d<dim; ++d) { - for (ibset::const_iterator - bi = notrefined.begin(); bi != notrefined.end(); ++bi) - { - if (hh.refcent == vertex_centered) { - enlarged[d] |= (*bi).expand(ivect::dir(d), ivect::dir(d)); - } else { - enlarged[d] |= *bi; - } - } - enlarged[d].normalize(); - } - - // Intersect with the union of refined regions - ibset boundaries[dim]; - for (int d=0; d<dim; ++d) { - boundaries[d] = active & enlarged[d]; - boundaries[d].normalize(); - } - - // Subtract the boundaries from the refined region - ibset refined = active; - for (int d=0; d<dim; ++d) { - refined -= boundaries[d]; - } - refined.normalize(); - - - // Set prolongation boundaries of this level { BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { DECLARE_CCTK_ARGUMENTS; - ibbox const & ext - = dd.boxes.at(mglevel).at(reflevel).at(component).exterior; + ibbox const & ext = + dd.light_boxes.AT(mglevel).AT(reflevel).AT(component).exterior; + ibset const & active = + dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).active; + + vect<ibset,dim> const & boundaries = + dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).prolongation_boundaries; - ibset notactive = ext - active; - notactive.normalize(); + ibset const notactive = ext - active; for (int d=0; d<dim; ++d) { assert ((notactive & boundaries[d]).empty()); } - for (ibset::const_iterator bi = notactive.begin(); - bi != notactive.end(); - ++bi) - { - ibbox const& box = *bi; - assert (box <= ext); - assert (not box.empty()); - - ivect imin, imax; - ibbox2iminimax (ext, box, ivect::ref(cctk_lsh), imin, imax); + LOOP_OVER_BSET (cctkGH, notactive, box, imin, imax) { if (verbose) { ostringstream buf; @@ -175,41 +94,31 @@ namespace CarpetMask { weight[ind] = 0.0; } LC_ENDLOOP3(CarpetMaskSetup_buffers); - } // for box + } END_LOOP_OVER_BSET; for (int d=0; d<dim; ++d) { - for (ibset::const_iterator bi = boundaries[d].begin(); - bi != boundaries[d].end(); - ++bi) - { + LOOP_OVER_BSET (cctkGH, boundaries[d], box, imin, imax) { - ibbox const & box = (*bi) & ext; - if (not box.empty()) { - - ivect imin, imax; - ibbox2iminimax (ext, box, ivect::ref(cctk_lsh), imin, imax); - - if (verbose) { - ostringstream buf; - buf << "Setting prolongation boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione; - CCTK_INFO (buf.str().c_str()); - } - - // Set weight on the boundary to 1/2 - assert (dim == 3); + if (verbose) { + ostringstream buf; + buf << "Setting prolongation boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione; + CCTK_INFO (buf.str().c_str()); + } + + // Set weight on the boundary to 1/2 + assert (dim == 3); #pragma omp parallel - LC_LOOP3(CarpetMaskSetup_prolongation, - i,j,k, - imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], - cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) - { - int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - weight[ind] *= 0.5; - } LC_ENDLOOP3(CarpetMaskSetup_prolongation); - - } // if box not empty + LC_LOOP3(CarpetMaskSetup_prolongation, + i,j,k, + imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + weight[ind] *= 0.5; + } LC_ENDLOOP3(CarpetMaskSetup_prolongation); + - } // for box + } END_LOOP_OVER_BSET; } // for d } END_LOCAL_COMPONENT_LOOP; @@ -231,41 +140,32 @@ namespace CarpetMask { DECLARE_CCTK_ARGUMENTS; - ibbox const & ext - = dd.boxes.at(mglevel).at(reflevel).at(component).exterior; + ibset const & refined = + dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).restricted_region; + vect<ibset,dim> const & boundaries = + dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).restriction_boundaries; - for (ibset::const_iterator bi = refined.begin(); - bi != refined.end(); - ++bi) - { + LOOP_OVER_BSET (cctkGH, refined, box, imin, imax) { - ibbox const & box = (*bi).contracted_for(ext) & ext; - if (not box.empty()) { - - ivect imin, imax; - ibbox2iminimax (ext, box, ivect::ref(cctk_lsh), imin, imax); - - if (verbose) { - ostringstream buf; - buf << "Setting restricted region on level " << reflevel << " to weight 0: " << imin << ":" << imax-ione; - CCTK_INFO (buf.str().c_str()); - } - - // Set weight in the restricted region to 0 - assert (dim == 3); + if (verbose) { + ostringstream buf; + buf << "Setting restricted region on level " << reflevel << " to weight 0: " << imin << ":" << imax-ione; + CCTK_INFO (buf.str().c_str()); + } + + // Set weight in the restricted region to 0 + assert (dim == 3); #pragma omp parallel - LC_LOOP3(CarpetMaskSetup_restriction, - i,j,k, - imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], - cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) - { - int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - weight[ind] = 0; - } LC_ENDLOOP3(CarpetMaskSetup_restriction); - - } // if box not empty - - } // for box + LC_LOOP3(CarpetMaskSetup_restriction, + i,j,k, + imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + weight[ind] = 0; + } LC_ENDLOOP3(CarpetMaskSetup_restriction); + + } END_LOOP_OVER_BSET; assert (dim == 3); vector<int> mask (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]); @@ -282,41 +182,30 @@ namespace CarpetMask { } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_init); for (int d=0; d<dim; ++d) { - for (ibset::const_iterator bi = boundaries[d].begin(); - bi != boundaries[d].end(); - ++bi) - { + LOOP_OVER_BSET (cctkGH, boundaries[d], box, imin, imax) { - ibbox const & box = (*bi).contracted_for(ext) & ext; - if (not box.empty()) { - - ivect imin, imax; - ibbox2iminimax (ext, box, ivect::ref(cctk_lsh), imin, imax); - - if (verbose) { - ostringstream buf; - buf << "Setting restriction boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione; - CCTK_INFO (buf.str().c_str()); - } - - // Set weight on the boundary to 1/2 - assert (dim == 3); + if (verbose) { + ostringstream buf; + buf << "Setting restriction boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione; + CCTK_INFO (buf.str().c_str()); + } + + // Set weight on the boundary to 1/2 + assert (dim == 3); #pragma omp parallel - LC_LOOP3(CarpetMaskSetup_restriction_boundary_partial, - i,j,k, - imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], - cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) - { - int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - if (mask[ind] == 0) { - mask[ind] = 1; - } - mask[ind] *= 2; - } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_partial); - - } // if box not empty + LC_LOOP3(CarpetMaskSetup_restriction_boundary_partial, + i,j,k, + imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + if (mask[ind] == 0) { + mask[ind] = 1; + } + mask[ind] *= 2; + } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_partial); - } // for box + } END_LOOP_OVER_BSET; } // for d assert (dim == 3); |