[Users] Carpet: using baseextents?

[Scott Hawley]

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 at cct.lsu.edu> wrote:

>On 2013-09-12, at 20:57 , Scott Hawley <scott.hawley at 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
>
>-- 
>Erik Schnetter <schnetter at 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/.
>