Thanks Erik. Glad I asked.
Since I *do* have multiple patches, how do I get the "map" patch number?
(Just referencing "map" didn't work, and checking BEGIN_MAP_LOOP in Carpet/src/modes.hh didn't enlighten me.)
-Scott
On 9/12/13 8:17 PM, "Erik Schnetter" schnetter@cct.lsu.edu wrote:
On 2013-09-12, at 20:57 , Scott Hawley scott.hawley@belmont.edu wrote:
Goal: I'm trying to determine if a given (fine) grid vertex coincides with a coarser grid point.
Frank pointed me to baseextents, and I've come up with the following snippet of code... Just posting this for a sanity check: Am I on the right track? Is there a more robust way one should do this?
My code relies on "baseextents" being the "full" extent of a given refinement level, apart from any domain decomposition that occurs when running with multiple MPI processors. (Well, and it relies on the fact that fine grid boundaries must lie along coarse grid points.)
The parenthesized remark doesn't sound quite right. See below, coarse_ilo.
BEGIN_MAP_LOOP (cctkGH, CCTK_GF) { BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
DECLARE_CCTK_ARGUMENTS; assert (maps == 1); int const m = 0; // <--- Really? m is never > 0?You are using "m" for two things. One of these should rather be called "ml", a "multigrid-level" that is currently unused, and is thus always 0. The other is the "map", i.e. the patch number; this is 0 unless you have multiple patches.
assert (mglevels == 1); int const rl = reflevel; ivect const ilo =vhh.at(m)->baseextents.at(m).at(rl).lower();
The second .at(m) should be .at(ml) instead.
ivect const istr =vhh.at(m)->baseextents.at(m).at(rl).stride();
Same here.
int ref_ratio = 2; // <-- Do we support refinementratios > 2?
Nothing uses this at the moment; it is untested and thus likely not working.
ivect coarse_str = istr * ref_ratio;Better: vhh.at(ml)->baseextents.at(m).at(rl-1).stride();
Also: coarse_ilo = vhh.at(m)->baseextents.at(ml).at(rl-1).lower();
for (CCTK_INT k = 0 ; k < cctk_lsh[2] ; k=k+1) { for (CCTK_INT j = 0 ; j < cctk_lsh[1] ; j=j+1) { for (CCTK_INT i = 0 ; i < cctk_lsh[0] ; i=i+1) { // does this vertex coincide with a coarser gridvertex? if ( ( 0 == (i - ilo[0]) % coarse_str[0]) && ( 0 == (j - ilo[1]) % coarse_str[1]) && ( 0 == (k - ilo[2]) % coarse_str[2])) {
This doesn't seem quite right. ilo[0] is not in the same "units" as i; instead, it is in the same "units" as str.
Alternatively: Check whether the equation ilo[0] + i * str[0] == coarse_ilo[0] + coarse_i * coarse_str[0] can be solved for coarse_i, i.e. check whether
(ilo[0] + i * str[0] - coarse_ilo[0]) % coarse_str[0] == 0
and similarly for the other two directions.
CCTK_VInfo(CCTK_THORNSTRING," We areover a CG vertex!"); }
} // end loop over i } // end loop over j } // end loop over k } END_LOCAL_COMPONENT_LOOP;} END_MAP_LOOP;
Thank you very much,
To test this, you can output the coordinate of the current grid point, and then check (manually) whether this coordinate exists on the next coarser grid.
-erik
-- Erik Schnetter schnetter@cct.lsu.edu http://www.perimeterinstitute.ca/personal/eschnetter/
My email is as private as my paper mail. I therefore support encrypting and signing email messages. Get my PGP key from http://pgp.mit.edu/.