Re: [Users] Carpet: using baseextents?

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 refinement ratios > 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 grid vertex?
                  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 are over

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