Scott
I notice that you access the variable baseextents. Instead, consider calling the function "baseextent", as in baseextent(ml,rl).
To get the extent of a particular component, write
extent(ml, rl, c)
where "c" is the component number. The "local" grid has component number "Carpet::component". This would be
extent(mglevel, reflevel, component)
which returns a bbox where lower(), upper(), and stride() define the current component.
-erik
On 2013-09-13, at 0:12 , Scott Hawley scott.hawley@belmont.edu wrote:
Erik,
I notice there's not "AT(c)" in the lines for the lower() calls, which would tell it which component grid I want the extents of.
The result is that both of those lower() lines return [0,0,0]. This is for the finest levels of multiply-nested grids. Similarly, if I replace "lower" with "upper", I get [768,768,768] -- ie it returns the bounds for the entire domain, not the grid at that level.
How can we get the integer-coordinate values for the boundaries of the local grid?
The stride() lines work just fine, by the way.
Thanks, Scott
On Sep 12, 2013, at 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 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 overa 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/.