On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkel hfinkel@anl.gov wrote:

Could I also decrease the block size? I currently have CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? Is there a restriction based on the number of ghost points?

Yes, you can reduce the block size. I assume that both the regridding operation and the time evolution will become slower if you do that, because more blocks will have to be handled. This will also require more memory, since each block carries its own ghost zone overhead.

No, there is no such restriction.

-erik

On Wed, 2011-08-31 at 13:53 -0500, Hal Finkel wrote:

On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote:

...

On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkel hfinkel@anl.gov wrote:

...

Could I also decrease the block size? I currently have CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? Is there a restriction based on the number of ghost points?

Yes, you can reduce the block size. I assume that both the regridding operation and the time evolution will become slower if you do that, because more blocks will have to be handled. This will also require more memory, since each block carries its own ghost zone overhead.

When I try reducing the block size from 4 to 3, it exits with:

build/CarpetLib/th.cc:107: void th::regrid(): Assertion `times.at(ml).at(rl).at(tl) < times.at(ml).at(rl).at(tl-1)' failed. p0_31544: p4_error: interrupt SIGx: 6

This happens with a larger box as well. It seems to have something to do with setting: Carpet::refine_timestep = yes.

-Hal

A block size of 2 is either really slow or it hangs.

I'll try a larger box for now.

Thanks again, Hal

...

No, there is no such restriction.

-erik

On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkel hfinkel@anl.gov wrote:

On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote:

...

On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkel hfinkel@anl.gov wrote:

...

Could I also decrease the block size? I currently have CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? Is there a restriction based on the number of ghost points?

Yes, you can reduce the block size. I assume that both the regridding operation and the time evolution will become slower if you do that, because more blocks will have to be handled.

Regardless of what I do, once we get past the first coarse time step, the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] Regridding map 0...".

Overall, it is in dh::regrid(do_init=true). It spends most of its time in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ nsi). The normalize() function does exit, however, so it is not hanging in that function.

The core problem seems to be that it takes a long time to execute: boxes  = boxes .shift(-dir) - boxes; in dh::regrid(do_init=true). Probably because boxes has 129064 elements. The coarse grid is now only 30^3 and I've left the regrid box size at 4. I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 ~ 420 refinement regions.

What is the best way to figure out what is going on?

Hal

Yes, this function is very slow. I did not expect it to be prohibitively slow. Are you compiling with optimisation enabled?

The bboxset represents the set of refined regions, and it is internally represented as a list of bboxes (regions). Carpet performs set operations on these (intersection, union, complement, etc.) to determine the communication schedule, i.e. which ghost zones of which bbox need to be filled from which other bbox. Unfortunately, the algorithm used for this is O(n^2) in the number of refined regions, and set operations when implemented via lists themselves are O(n^2) in the set size, leading to a rather unfortunate overall complexity. The only cure is to reduce the number of bboxes (make them larger) and to regrid fewer times.

We are running on several hundred MPI processes, and are therefore dealing with set sizes of several hundred in production mode. (Carpet's domain decomposition splits boxes, so that each process receives a set of boxes.) It may be that the irregularity of the AMR box distribution leads to a worse case than our domain decomposition.

We have implemented / are implementing set operations based on a tree datastructure, which should be much more efficient. This is the work of Ashley Zebrowski, a graduate student at LSU who just tried (unsuccessfully -- a hurricane was in the way) to attend the ParCo 2011 conference in Belgium to present this work.

Frank, could you give a brief update of the status of Ashley's work?

-erik

On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkel hfinkel@anl.gov wrote:

On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkel hfinkel@anl.gov wrote:

...

On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote:

...

On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkel hfinkel@anl.gov wrote:

...

Could I also decrease the block size? I currently have CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? Is there a restriction based on the number of ghost points?

Yes, you can reduce the block size. I assume that both the regridding operation and the time evolution will become slower if you do that, because more blocks will have to be handled.

Regardless of what I do, once we get past the first coarse time step, the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] Regridding map 0...".

Overall, it is in dh::regrid(do_init=true). It spends most of its time in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ nsi). The normalize() function does exit, however, so it is not hanging in that function.

The core problem seems to be that it takes a long time to execute: boxes  = boxes .shift(-dir) - boxes; in dh::regrid(do_init=true). Probably because boxes has 129064 elements. The coarse grid is now only 30^3 and I've left the regrid box size at 4. I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 ~ 420 refinement regions.

What is the best way to figure out what is going on?

Hal

Yes, this function is very slow. I did not expect it to be prohibitively slow. Are you compiling with optimisation enabled?

I've tried with optimizations enabled (and without for debugging).

...

The bboxset represents the set of refined regions, and it is internally represented as a list of bboxes (regions). Carpet performs set operations on these (intersection, union, complement, etc.) to determine the communication schedule, i.e. which ghost zones of which bbox need to be filled from which other bbox. Unfortunately, the algorithm used for this is O(n^2) in the number of refined regions, and set operations when implemented via lists themselves are O(n^2) in the set size, leading to a rather unfortunate overall complexity. The only cure is to reduce the number of bboxes (make them larger) and to regrid fewer times.

This is what I suspected, but nevertheless, is there something wrong? How many boxes do you expect that I should have? The reason that it does not finish, even with optimizations, is that there are 129K boxes in the loop (that's at least 16 billion box normalizations?).

The coarse grid is only 30^3, and the regrid box size is 4, so at maximum, there should be ~400 level one boxes. Even if some of those have level 2 boxes, I don't understand how there could be 129K boxes.

The refinement structure itself should have one bbox per refined 4^3 box, and both CarpetRegrid2 and CarpetLib would try to combine these into fewer boxes where possible, i.e. where one can form rectangles or larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on level one.

Internally, one of the expensive operations is to determine which points are buffer points. But before I explain this algorithm one important question: Are you using buffer points? Those become expensive, and you should be able to avoid them. First, you should be able to just switch off this parameter; if this leads to loss of convergence, you should be able to rewrite your equations into a first-order system (only first spatial derivatives by adding another evolved variable), which will most likely be stable. Not using buffer zones should make things much cheaper.

-erik

On Thu, Sep 1, 2011 at 3:05 PM, Hal Finkel hfinkel@anl.gov wrote:

On Thu, 2011-09-01 at 14:25 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkel hfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkel hfinkel@anl.gov wrote:

...

On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote:

...

On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkel hfinkel@anl.gov wrote: > Could I also decrease the block size? I currently have > CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? > Is there a restriction based on the number of ghost points?

Yes, you can reduce the block size. I assume that both the regridding operation and the time evolution will become slower if you do that, because more blocks will have to be handled.

Regardless of what I do, once we get past the first coarse time step, the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] Regridding map 0...".

Overall, it is in dh::regrid(do_init=true). It spends most of its time in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ nsi). The normalize() function does exit, however, so it is not hanging in that function.

The core problem seems to be that it takes a long time to execute: boxes  = boxes .shift(-dir) - boxes; in dh::regrid(do_init=true). Probably because boxes has 129064 elements. The coarse grid is now only 30^3 and I've left the regrid box size at 4. I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 ~ 420 refinement regions.

What is the best way to figure out what is going on?

Hal

Yes, this function is very slow. I did not expect it to be prohibitively slow. Are you compiling with optimisation enabled?

I've tried with optimizations enabled (and without for debugging).

...

The bboxset represents the set of refined regions, and it is internally represented as a list of bboxes (regions). Carpet performs set operations on these (intersection, union, complement, etc.) to determine the communication schedule, i.e. which ghost zones of which bbox need to be filled from which other bbox. Unfortunately, the algorithm used for this is O(n^2) in the number of refined regions, and set operations when implemented via lists themselves are O(n^2) in the set size, leading to a rather unfortunate overall complexity. The only cure is to reduce the number of bboxes (make them larger) and to regrid fewer times.

This is what I suspected, but nevertheless, is there something wrong? How many boxes do you expect that I should have? The reason that it does not finish, even with optimizations, is that there are 129K boxes in the loop (that's at least 16 billion box normalizations?).

The coarse grid is only 30^3, and the regrid box size is 4, so at maximum, there should be ~400 level one boxes. Even if some of those have level 2 boxes, I don't understand how there could be 129K boxes.

The refinement structure itself should have one bbox per refined 4^3 box, and both CarpetRegrid2 and CarpetLib would try to combine these into fewer boxes where possible, i.e. where one can form rectangles or larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on level one.

That makes sense. I think that there is a bug somewhere which is causing the box set to be much too big. Furthermore, it does not happen on every run, only sometimes. When it does not happen, I hit another bug after a few coarse timesteps:

I get a range-check exception from std::vector in a call to: gh::get_local_component (rl=1, c=8) the problem is that this returns: local_components_.AT(rl).AT(c); and local_components_[1].size() is 8 The call to get_local_component is coming from ggf::transfer_from_all at: int const lc2 = h.get_local_component(rl2,c2); where c2 is from psend.component. So it looks like there is an off-by-one error somewhere.

Very strange. This code should be quite solid by now. psend is set in the file dh.cc in thorn Carpet/CarpetLib; there is one (large) routine that calculates the communication schedule. Some of the indexing errors there in the past included confusing the number of components on different refinement levels, which led to indexing errors such as the one you describe.

...

Internally, one of the expensive operations is to determine which points are buffer points. But before I explain this algorithm one important question: Are you using buffer points? Those become expensive, and you should be able to avoid them. First, you should be able to just switch off this parameter; if this leads to loss of convergence, you should be able to rewrite your equations into a first-order system (only first spatial derivatives by adding another evolved variable), which will most likely be stable. Not using buffer zones should make things much cheaper.

I am computing second spatial derivatives, and while I am curious about the algorithm, etc. I am not sure this has anything to do with the problem. I've changed the threshold so that I'm not generating level_mask values >= 2 and I still sometimes end up in that loop with 170K boxes. That is clearly just wrong. If you can tell me where the boxes are generated (in which functions), I can add some printf statements and breakpoints and I'll try and see where they are all coming from.

Buffer zones increase the refined region. Try Carpet::use_buffer_zones=no, which is also the default.

The level mask is interpreted in the file amr.cc, function evaluate_level_mask in thorn Carpet/CarpetRegrid2. These boxes are then post-processed in file regrid.cc, lines 315 to 415. The individual post-processing steps are defined in file property.cc. The post-processing enforces certain grid structure properties, such as e.g. proper nesting, certain symmetries (this is also where periodicity would be enforced), or adding buffer zones.

-erik

Thanks again, Hal

...

-erik

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

On Thu, Sep 1, 2011 at 3:56 PM, Hal Finkel hfinkel@anl.gov wrote:

On Thu, 2011-09-01 at 15:16 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 3:05 PM, Hal Finkel hfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 14:25 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkel hfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkel hfinkel@anl.gov wrote: > On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote: >> On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkel hfinkel@anl.gov wrote: >> > Could I also decrease the block size? I currently have >> > CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? >> > Is there a restriction based on the number of ghost points? >> >> Yes, you can reduce the block size. I assume that both the regridding >> operation and the time evolution will become slower if you do that, >> because more blocks will have to be handled. > > Regardless of what I do, once we get past the first coarse time step, > the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] > Regridding map 0...". > > Overall, it is in dh::regrid(do_init=true). It spends most of its time > in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: > for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ > nsi). The normalize() function does exit, however, so it is not hanging > in that function. > > The core problem seems to be that it takes a long time to execute: > boxes  = boxes .shift(-dir) - boxes; > in dh::regrid(do_init=true). Probably because boxes has 129064 elements. > The coarse grid is now only 30^3 and I've left the regrid box size at 4. > I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 > ~ 420 refinement regions. > > What is the best way to figure out what is going on?

Hal

Yes, this function is very slow. I did not expect it to be prohibitively slow. Are you compiling with optimisation enabled?

I've tried with optimizations enabled (and without for debugging).

...

The bboxset represents the set of refined regions, and it is internally represented as a list of bboxes (regions). Carpet performs set operations on these (intersection, union, complement, etc.) to determine the communication schedule, i.e. which ghost zones of which bbox need to be filled from which other bbox. Unfortunately, the algorithm used for this is O(n^2) in the number of refined regions, and set operations when implemented via lists themselves are O(n^2) in the set size, leading to a rather unfortunate overall complexity. The only cure is to reduce the number of bboxes (make them larger) and to regrid fewer times.

This is what I suspected, but nevertheless, is there something wrong? How many boxes do you expect that I should have? The reason that it does not finish, even with optimizations, is that there are 129K boxes in the loop (that's at least 16 billion box normalizations?).

The coarse grid is only 30^3, and the regrid box size is 4, so at maximum, there should be ~400 level one boxes. Even if some of those have level 2 boxes, I don't understand how there could be 129K boxes.

The refinement structure itself should have one bbox per refined 4^3 box, and both CarpetRegrid2 and CarpetLib would try to combine these into fewer boxes where possible, i.e. where one can form rectangles or larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on level one.

That makes sense. I think that there is a bug somewhere which is causing the box set to be much too big. Furthermore, it does not happen on every run, only sometimes. When it does not happen, I hit another bug after a few coarse timesteps:

I get a range-check exception from std::vector in a call to: gh::get_local_component (rl=1, c=8) the problem is that this returns: local_components_.AT(rl).AT(c); and local_components_[1].size() is 8 The call to get_local_component is coming from ggf::transfer_from_all at: int const lc2 = h.get_local_component(rl2,c2); where c2 is from psend.component. So it looks like there is an off-by-one error somewhere.

Very strange. This code should be quite solid by now. psend is set in the file dh.cc in thorn Carpet/CarpetLib; there is one (large) routine that calculates the communication schedule. Some of the indexing errors there in the past included confusing the number of components on different refinement levels, which led to indexing errors such as the one you describe.

...

...

Internally, one of the expensive operations is to determine which points are buffer points. But before I explain this algorithm one important question: Are you using buffer points? Those become expensive, and you should be able to avoid them. First, you should be able to just switch off this parameter; if this leads to loss of convergence, you should be able to rewrite your equations into a first-order system (only first spatial derivatives by adding another evolved variable), which will most likely be stable. Not using buffer zones should make things much cheaper.

I am computing second spatial derivatives, and while I am curious about the algorithm, etc. I am not sure this has anything to do with the problem. I've changed the threshold so that I'm not generating level_mask values >= 2 and I still sometimes end up in that loop with 170K boxes. That is clearly just wrong. If you can tell me where the boxes are generated (in which functions), I can add some printf statements and breakpoints and I'll try and see where they are all coming from.

Buffer zones increase the refined region. Try Carpet::use_buffer_zones=no, which is also the default.

I did not have that turned on.

What is the different between boxes.size() and boxes.setsize()? Perhaps I was misinterpreting these: (gdb) p boxes.setsize() $1 = 24 (gdb) p boxes.size() $2 = 175175 (gdb)

The set size is the number of bboxes in the list, the size is the number of grid points in all bboxes combined. The algorithm cost depends on the setsize only. Good; this clears up part of the problem.

So maybe I don't have too many boxes, maybe it is just hanging sometimes (and then otherwise crashing with that indexing error).

In terms of non-default Carpet parameters, I did have: regrid_in_level_mode = no (to suppress some warning early on). Changing that back to its default value of "yes" does not seem to make a difference.

This parameter only makes a difference if you have multiple maps, which you don't. I would nevertheless leave it on, since that is the default, and may traverse code paths that are more well tested. Which warning do you see? The warning enabled by regrid_during_initialisation should be harmless.

I now have: Carpet::verbose           = yes Carpet::veryverbose       = yes Carpet::schedule_barriers = no Carpet::storage_verbose   = no CarpetLib::output_bboxes  = no

You can use this parameter to output the gory details of the grid structure. I use this for debugging, but it really contains a lot of detail.

Carpet::output_timers_every      = 512 CarpetLib::print_timestats_every = 512 CarpetLib::print_memstats_every  = 512 Carpet::domain_from_coordbase = yes Carpet::prolongation_order_time =  2 Carpet::prolongation_order_space = 3 driver::ghost_size = 3 Carpet::poison_new_timelevels = yes CarpetLib::poison_new_memory  = yes Carpet::init_3_timelevels = yes

I would leave this off; Carpet::init_each_timelevel or Carpet::init_fill_timelevels would be better. init_3_timelevels performs some time evolution steps forwards and backwards to fill the past timelevels; before we debugged the remainder, I would be afraid of regridding too often. init_each_timelevel calls the initial data routine multiple times with different cctk_time, init_fill_timelevels copies the current timelevel to the past timelevels.

Carpet::max_refinement_levels = 10 Carpet::refine_timestep = no Carpet::use_buffer_zones = no CarpetRegrid2::verbose                 = yes CarpetRegrid2::veryverbose             = yes CarpetRegrid2::regrid_every = 512 CarpetRegrid2::adaptive_refinement = yes CarpetRegrid2::adaptive_block_size = 4 CarpetRegrid2::add_levels_automatically = yes

These are all good.

-erik

Thanks again, Hal

...

The level mask is interpreted in the file amr.cc, function evaluate_level_mask in thorn Carpet/CarpetRegrid2. These boxes are then post-processed in file regrid.cc, lines 315 to 415. The individual post-processing steps are defined in file property.cc. The post-processing enforces certain grid structure properties, such as e.g. proper nesting, certain symmetries (this is also where periodicity would be enforced), or adding buffer zones.

-erik

...

Thanks again, Hal

...

-erik

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

On Thu, Sep 1, 2011 at 4:29 PM, Hal Finkel hfinkel@anl.gov wrote:

On Thu, 2011-09-01 at 16:14 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 3:56 PM, Hal Finkel hfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 15:16 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 3:05 PM, Hal Finkel hfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 14:25 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkel hfinkel@anl.gov wrote: > On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote: >> On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkel hfinkel@anl.gov wrote: >> > On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote: >> >> On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkel hfinkel@anl.gov wrote: >> >> > Could I also decrease the block size? I currently have >> >> > CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? >> >> > Is there a restriction based on the number of ghost points? >> >> >> >> Yes, you can reduce the block size. I assume that both the regridding >> >> operation and the time evolution will become slower if you do that, >> >> because more blocks will have to be handled. >> > >> > Regardless of what I do, once we get past the first coarse time step, >> > the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] >> > Regridding map 0...". >> > >> > Overall, it is in dh::regrid(do_init=true). It spends most of its time >> > in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: >> > for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ >> > nsi). The normalize() function does exit, however, so it is not hanging >> > in that function. >> > >> > The core problem seems to be that it takes a long time to execute: >> > boxes  = boxes .shift(-dir) - boxes; >> > in dh::regrid(do_init=true). Probably because boxes has 129064 elements. >> > The coarse grid is now only 30^3 and I've left the regrid box size at 4. >> > I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 >> > ~ 420 refinement regions. >> > >> > What is the best way to figure out what is going on? >> >> Hal >> >> Yes, this function is very slow. I did not expect it to be >> prohibitively slow. Are you compiling with optimisation enabled? > > I've tried with optimizations enabled (and without for debugging). > >> >> The bboxset represents the set of refined regions, and it is >> internally represented as a list of bboxes (regions). Carpet performs >> set operations on these (intersection, union, complement, etc.) to >> determine the communication schedule, i.e. which ghost zones of which >> bbox need to be filled from which other bbox. Unfortunately, the >> algorithm used for this is O(n^2) in the number of refined regions, >> and set operations when implemented via lists themselves are O(n^2) in >> the set size, leading to a rather unfortunate overall complexity. The >> only cure is to reduce the number of bboxes (make them larger) and to >> regrid fewer times. > > This is what I suspected, but nevertheless, is there something wrong? > How many boxes do you expect that I should have? The reason that it does > not finish, even with optimizations, is that there are 129K boxes in the > loop (that's at least 16 billion box normalizations?). > > The coarse grid is only 30^3, and the regrid box size is 4, so at > maximum, there should be ~400 level one boxes. Even if some of those > have level 2 boxes, I don't understand how there could be 129K boxes.

The refinement structure itself should have one bbox per refined 4^3 box, and both CarpetRegrid2 and CarpetLib would try to combine these into fewer boxes where possible, i.e. where one can form rectangles or larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on level one.

That makes sense. I think that there is a bug somewhere which is causing the box set to be much too big. Furthermore, it does not happen on every run, only sometimes. When it does not happen, I hit another bug after a few coarse timesteps:

I get a range-check exception from std::vector in a call to: gh::get_local_component (rl=1, c=8) the problem is that this returns: local_components_.AT(rl).AT(c); and local_components_[1].size() is 8 The call to get_local_component is coming from ggf::transfer_from_all at: int const lc2 = h.get_local_component(rl2,c2); where c2 is from psend.component. So it looks like there is an off-by-one error somewhere.

Very strange. This code should be quite solid by now. psend is set in the file dh.cc in thorn Carpet/CarpetLib; there is one (large) routine that calculates the communication schedule. Some of the indexing errors there in the past included confusing the number of components on different refinement levels, which led to indexing errors such as the one you describe.

...

...

Internally, one of the expensive operations is to determine which points are buffer points. But before I explain this algorithm one important question: Are you using buffer points? Those become expensive, and you should be able to avoid them. First, you should be able to just switch off this parameter; if this leads to loss of convergence, you should be able to rewrite your equations into a first-order system (only first spatial derivatives by adding another evolved variable), which will most likely be stable. Not using buffer zones should make things much cheaper.

I am computing second spatial derivatives, and while I am curious about the algorithm, etc. I am not sure this has anything to do with the problem. I've changed the threshold so that I'm not generating level_mask values >= 2 and I still sometimes end up in that loop with 170K boxes. That is clearly just wrong. If you can tell me where the boxes are generated (in which functions), I can add some printf statements and breakpoints and I'll try and see where they are all coming from.

Buffer zones increase the refined region. Try Carpet::use_buffer_zones=no, which is also the default.

I did not have that turned on.

What is the different between boxes.size() and boxes.setsize()? Perhaps I was misinterpreting these: (gdb) p boxes.setsize() $1 = 24 (gdb) p boxes.size() $2 = 175175 (gdb)

The set size is the number of bboxes in the list, the size is the number of grid points in all bboxes combined. The algorithm cost depends on the setsize only. Good; this clears up part of the problem.

...

So maybe I don't have too many boxes, maybe it is just hanging sometimes (and then otherwise crashing with that indexing error).

In terms of non-default Carpet parameters, I did have: regrid_in_level_mode = no (to suppress some warning early on). Changing that back to its default value of "yes" does not seem to make a difference.

This parameter only makes a difference if you have multiple maps, which you don't. I would nevertheless leave it on, since that is the default, and may traverse code paths that are more well tested. Which warning do you see? The warning enabled by regrid_during_initialisation should be harmless.

That was the one. I'll ignore it if I turn that back on.

...

...

I now have: Carpet::verbose           = yes Carpet::veryverbose       = yes Carpet::schedule_barriers = no Carpet::storage_verbose   = no CarpetLib::output_bboxes  = no

You can use this parameter to output the gory details of the grid structure. I use this for debugging, but it really contains a lot of detail.

...

Carpet::output_timers_every      = 512 CarpetLib::print_timestats_every = 512 CarpetLib::print_memstats_every  = 512 Carpet::domain_from_coordbase = yes Carpet::prolongation_order_time =  2 Carpet::prolongation_order_space = 3 driver::ghost_size = 3 Carpet::poison_new_timelevels = yes CarpetLib::poison_new_memory  = yes Carpet::init_3_timelevels = yes

I would leave this off; Carpet::init_each_timelevel or Carpet::init_fill_timelevels would be better. init_3_timelevels performs some time evolution steps forwards and backwards to fill the past timelevels; before we debugged the remainder, I would be afraid of regridding too often. init_each_timelevel calls the initial data routine multiple times with different cctk_time, init_fill_timelevels copies the current timelevel to the past timelevels.

Does this have anything to do with the timelevels=3 in my interface file? (This is there because the example from which I copied had Timelevels -> 3 in the Kranc driver code). Should I change it to 2?

For second order time interpolation (which you request), you need at least 3 timelevels. init_3_timelevels works only with 3 timelevels; init_each and init_fill work with any number of timelevels. You should keep 3 timelevels.

Thanks again, Hal

...

...

Carpet::max_refinement_levels = 10 Carpet::refine_timestep = no Carpet::use_buffer_zones = no CarpetRegrid2::verbose                 = yes CarpetRegrid2::veryverbose             = yes CarpetRegrid2::regrid_every = 512 CarpetRegrid2::adaptive_refinement = yes CarpetRegrid2::adaptive_block_size = 4 CarpetRegrid2::add_levels_automatically = yes

These are all good.

-erik

...

Thanks again, Hal

...

The level mask is interpreted in the file amr.cc, function evaluate_level_mask in thorn Carpet/CarpetRegrid2. These boxes are then post-processed in file regrid.cc, lines 315 to 415. The individual post-processing steps are defined in file property.cc. The post-processing enforces certain grid structure properties, such as e.g. proper nesting, certain symmetries (this is also where periodicity would be enforced), or adding buffer zones.

-erik

...

Thanks again, Hal

...

-erik

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

On Thu, 2011-09-01 at 15:16 -0400, Erik Schnetter wrote:

On Thu, Sep 1, 2011 at 3:05 PM, Hal Finkel hfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 14:25 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkel hfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkel hfinkel@anl.gov wrote:

...

On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote: > On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkel hfinkel@anl.gov wrote: > > Could I also decrease the block size? I currently have > > CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? > > Is there a restriction based on the number of ghost points? > > Yes, you can reduce the block size. I assume that both the regridding > operation and the time evolution will become slower if you do that, > because more blocks will have to be handled.

Regardless of what I do, once we get past the first coarse time step, the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] Regridding map 0...".

Overall, it is in dh::regrid(do_init=true). It spends most of its time in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ nsi). The normalize() function does exit, however, so it is not hanging in that function.

The core problem seems to be that it takes a long time to execute: boxes = boxes .shift(-dir) - boxes; in dh::regrid(do_init=true). Probably because boxes has 129064 elements. The coarse grid is now only 30^3 and I've left the regrid box size at 4. I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 ~ 420 refinement regions.

What is the best way to figure out what is going on?

Hal

Yes, this function is very slow. I did not expect it to be prohibitively slow. Are you compiling with optimisation enabled?

I've tried with optimizations enabled (and without for debugging).

...

The bboxset represents the set of refined regions, and it is internally represented as a list of bboxes (regions). Carpet performs set operations on these (intersection, union, complement, etc.) to determine the communication schedule, i.e. which ghost zones of which bbox need to be filled from which other bbox. Unfortunately, the algorithm used for this is O(n^2) in the number of refined regions, and set operations when implemented via lists themselves are O(n^2) in the set size, leading to a rather unfortunate overall complexity. The only cure is to reduce the number of bboxes (make them larger) and to regrid fewer times.

This is what I suspected, but nevertheless, is there something wrong? How many boxes do you expect that I should have? The reason that it does not finish, even with optimizations, is that there are 129K boxes in the loop (that's at least 16 billion box normalizations?).

The coarse grid is only 30^3, and the regrid box size is 4, so at maximum, there should be ~400 level one boxes. Even if some of those have level 2 boxes, I don't understand how there could be 129K boxes.

The refinement structure itself should have one bbox per refined 4^3 box, and both CarpetRegrid2 and CarpetLib would try to combine these into fewer boxes where possible, i.e. where one can form rectangles or larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on level one.

That makes sense. I think that there is a bug somewhere which is causing the box set to be much too big. Furthermore, it does not happen on every run, only sometimes. When it does not happen, I hit another bug after a few coarse timesteps:

I get a range-check exception from std::vector in a call to: gh::get_local_component (rl=1, c=8) the problem is that this returns: local_components_.AT(rl).AT(c); and local_components_[1].size() is 8 The call to get_local_component is coming from ggf::transfer_from_all at: int const lc2 = h.get_local_component(rl2,c2); where c2 is from psend.component. So it looks like there is an off-by-one error somewhere.

Very strange. This code should be quite solid by now. psend is set in the file dh.cc in thorn Carpet/CarpetLib; there is one (large) routine that calculates the communication schedule. Some of the indexing errors there in the past included confusing the number of components on different refinement levels, which led to indexing errors such as the one you describe.

The bad component numbers are not coming from: preg.component = tmpncomps.AT(m)++; in Carpet/src/Recompose.cc

Where else are the component numbers assigned?

Thanks again, Hal

...

...

Internally, one of the expensive operations is to determine which points are buffer points. But before I explain this algorithm one important question: Are you using buffer points? Those become expensive, and you should be able to avoid them. First, you should be able to just switch off this parameter; if this leads to loss of convergence, you should be able to rewrite your equations into a first-order system (only first spatial derivatives by adding another evolved variable), which will most likely be stable. Not using buffer zones should make things much cheaper.

I am computing second spatial derivatives, and while I am curious about the algorithm, etc. I am not sure this has anything to do with the problem. I've changed the threshold so that I'm not generating level_mask values >= 2 and I still sometimes end up in that loop with 170K boxes. That is clearly just wrong. If you can tell me where the boxes are generated (in which functions), I can add some printf statements and breakpoints and I'll try and see where they are all coming from.

Buffer zones increase the refined region. Try Carpet::use_buffer_zones=no, which is also the default.

The level mask is interpreted in the file amr.cc, function evaluate_level_mask in thorn Carpet/CarpetRegrid2. These boxes are then post-processed in file regrid.cc, lines 315 to 415. The individual post-processing steps are defined in file property.cc. The post-processing enforces certain grid structure properties, such as e.g. proper nesting, certain symmetries (this is also where periodicity would be enforced), or adding buffer zones.

-erik

...

Thanks again, Hal

...

-erik

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

Erik, Hal,

Did you have any luck tracking down this error? I've just come back to this and am seeing the same error message; it appears to arise when two grids on a refined level merge, as in:

[3][0][0] exterior: [-0.010000,-0.020000,-0.020000] : [0.000000,0.020000,0.020000] : [0.001250,0.001250,0.001250] [3][0][1] exterior: [0.010000,-0.020000,-0.020000] : [0.040000,0.020000,0.020000] : [0.001250,0.001250,0.001250]

becomes

[3][0][0] exterior: [-0.010000,-0.020000,-0.020000] : [0.040000,0.020000,0.020000] : [0.001250,0.001250,0.001250]

It seems that the old grids are destroyed before the data is copied/populated in Recompose (either that or the old grid structure is not referred to in the data transfer).

Ian

On 07/09/11 01:03, Erik Schnetter wrote:

Hal

This is where numbers are assigned to components. The communication schedule decides which component needs to send data to which other component (which may be located on another process or not); this schedule is created for each refinement level independently, and may (if there is an error) refer to component numbers that don't exist. This schedule is set up in dh.cc.

Can you send me the example you are currently running (your source code and parameter file)? I will try to give it a try.

-erik

On Tue, Sep 6, 2011 at 7:44 PM, Hal Finkelhfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 15:16 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 3:05 PM, Hal Finkelhfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 14:25 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkelhfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote: > On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkelhfinkel@anl.gov wrote: >> On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote: >>> On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkelhfinkel@anl.gov wrote: >>>> Could I also decrease the block size? I currently have >>>> CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? >>>> Is there a restriction based on the number of ghost points? >>> Yes, you can reduce the block size. I assume that both the regridding >>> operation and the time evolution will become slower if you do that, >>> because more blocks will have to be handled. >> Regardless of what I do, once we get past the first coarse time step, >> the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] >> Regridding map 0...". >> >> Overall, it is in dh::regrid(do_init=true). It spends most of its time >> in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: >> for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ >> nsi). The normalize() function does exit, however, so it is not hanging >> in that function. >> >> The core problem seems to be that it takes a long time to execute: >> boxes = boxes .shift(-dir) - boxes; >> in dh::regrid(do_init=true). Probably because boxes has 129064 elements. >> The coarse grid is now only 30^3 and I've left the regrid box size at 4. >> I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 >> ~ 420 refinement regions. >> >> What is the best way to figure out what is going on? > Hal > > Yes, this function is very slow. I did not expect it to be > prohibitively slow. Are you compiling with optimisation enabled? I've tried with optimizations enabled (and without for debugging).

> The bboxset represents the set of refined regions, and it is > internally represented as a list of bboxes (regions). Carpet performs > set operations on these (intersection, union, complement, etc.) to > determine the communication schedule, i.e. which ghost zones of which > bbox need to be filled from which other bbox. Unfortunately, the > algorithm used for this is O(n^2) in the number of refined regions, > and set operations when implemented via lists themselves are O(n^2) in > the set size, leading to a rather unfortunate overall complexity. The > only cure is to reduce the number of bboxes (make them larger) and to > regrid fewer times. This is what I suspected, but nevertheless, is there something wrong? How many boxes do you expect that I should have? The reason that it does not finish, even with optimizations, is that there are 129K boxes in the loop (that's at least 16 billion box normalizations?).

The coarse grid is only 30^3, and the regrid box size is 4, so at maximum, there should be ~400 level one boxes. Even if some of those have level 2 boxes, I don't understand how there could be 129K boxes.

The refinement structure itself should have one bbox per refined 4^3 box, and both CarpetRegrid2 and CarpetLib would try to combine these into fewer boxes where possible, i.e. where one can form rectangles or larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on level one.

That makes sense. I think that there is a bug somewhere which is causing the box set to be much too big. Furthermore, it does not happen on every run, only sometimes. When it does not happen, I hit another bug after a few coarse timesteps:

I get a range-check exception from std::vector in a call to: gh::get_local_component (rl=1, c=8) the problem is that this returns: local_components_.AT(rl).AT(c); and local_components_[1].size() is 8 The call to get_local_component is coming from ggf::transfer_from_all at: int const lc2 = h.get_local_component(rl2,c2); where c2 is from psend.component. So it looks like there is an off-by-one error somewhere.

Very strange. This code should be quite solid by now. psend is set in the file dh.cc in thorn Carpet/CarpetLib; there is one (large) routine that calculates the communication schedule. Some of the indexing errors there in the past included confusing the number of components on different refinement levels, which led to indexing errors such as the one you describe.

The bad component numbers are not coming from: preg.component = tmpncomps.AT(m)++; in Carpet/src/Recompose.cc

Where else are the component numbers assigned?

Hal

Were you running with multiple processes?

-erik

On Tue, Sep 13, 2011 at 12:00 PM, Hal Finkel hfinkel@anl.gov wrote:

Ian,

I'd not found the problem yet. Hopefully this will be the hint we need to get this fixed.

Thanks again, Hal

On Tue, 2011-09-13 at 15:31 +0100, Ian Hawke wrote:

...

Erik, Hal,

Did you have any luck tracking down this error? I've just come back to this and am seeing the same error message; it appears to arise when two grids on a refined level merge, as in:

[3][0][0]   exterior: [-0.010000,-0.020000,-0.020000] : [0.000000,0.020000,0.020000] : [0.001250,0.001250,0.001250]     [3][0][1]   exterior: [0.010000,-0.020000,-0.020000] : [0.040000,0.020000,0.020000] : [0.001250,0.001250,0.001250]

becomes

[3][0][0]   exterior: [-0.010000,-0.020000,-0.020000] : [0.040000,0.020000,0.020000] : [0.001250,0.001250,0.001250]

It seems that the old grids are destroyed before the data is copied/populated in Recompose (either that or the old grid structure is not referred to in the data transfer).

Ian

On 07/09/11 01:03, Erik Schnetter wrote:

...

Hal

This is where numbers are assigned to components. The communication schedule decides which component needs to send data to which other component (which may be located on another process or not); this schedule is created for each refinement level independently, and may (if there is an error) refer to component numbers that don't exist. This schedule is set up in dh.cc.

Can you send me the example you are currently running (your source code and parameter file)? I will try to give it a try.

-erik

On Tue, Sep 6, 2011 at 7:44 PM, Hal Finkelhfinkel@anl.gov  wrote:

...

On Thu, 2011-09-01 at 15:16 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 3:05 PM, Hal Finkelhfinkel@anl.gov  wrote:

...

On Thu, 2011-09-01 at 14:25 -0400, Erik Schnetter wrote: > On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkelhfinkel@anl.gov  wrote: >> On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote: >>> On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkelhfinkel@anl.gov  wrote: >>>> On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote: >>>>> On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkelhfinkel@anl.gov  wrote: >>>>>> Could I also decrease the block size? I currently have >>>>>> CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? >>>>>> Is there a restriction based on the number of ghost points? >>>>> Yes, you can reduce the block size. I assume that both the regridding >>>>> operation and the time evolution will become slower if you do that, >>>>> because more blocks will have to be handled. >>>> Regardless of what I do, once we get past the first coarse time step, >>>> the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] >>>> Regridding map 0...". >>>> >>>> Overall, it is in dh::regrid(do_init=true). It spends most of its time >>>> in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: >>>> for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ >>>> nsi). The normalize() function does exit, however, so it is not hanging >>>> in that function. >>>> >>>> The core problem seems to be that it takes a long time to execute: >>>> boxes  = boxes .shift(-dir) - boxes; >>>> in dh::regrid(do_init=true). Probably because boxes has 129064 elements. >>>> The coarse grid is now only 30^3 and I've left the regrid box size at 4. >>>> I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 >>>> ~ 420 refinement regions. >>>> >>>> What is the best way to figure out what is going on? >>> Hal >>> >>> Yes, this function is very slow. I did not expect it to be >>> prohibitively slow. Are you compiling with optimisation enabled? >> I've tried with optimizations enabled (and without for debugging). >> >>> The bboxset represents the set of refined regions, and it is >>> internally represented as a list of bboxes (regions). Carpet performs >>> set operations on these (intersection, union, complement, etc.) to >>> determine the communication schedule, i.e. which ghost zones of which >>> bbox need to be filled from which other bbox. Unfortunately, the >>> algorithm used for this is O(n^2) in the number of refined regions, >>> and set operations when implemented via lists themselves are O(n^2) in >>> the set size, leading to a rather unfortunate overall complexity. The >>> only cure is to reduce the number of bboxes (make them larger) and to >>> regrid fewer times. >> This is what I suspected, but nevertheless, is there something wrong? >> How many boxes do you expect that I should have? The reason that it does >> not finish, even with optimizations, is that there are 129K boxes in the >> loop (that's at least 16 billion box normalizations?). >> >> The coarse grid is only 30^3, and the regrid box size is 4, so at >> maximum, there should be ~400 level one boxes. Even if some of those >> have level 2 boxes, I don't understand how there could be 129K boxes. > The refinement structure itself should have one bbox per refined 4^3 > box, and both CarpetRegrid2 and CarpetLib would try to combine these > into fewer boxes where possible, i.e. where one can form rectangles or > larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on > level one. That makes sense. I think that there is a bug somewhere which is causing the box set to be much too big. Furthermore, it does not happen on every run, only sometimes. When it does not happen, I hit another bug after a few coarse timesteps:

I get a range-check exception from std::vector in a call to: gh::get_local_component (rl=1, c=8) the problem is that this returns: local_components_.AT(rl).AT(c); and local_components_[1].size() is 8 The call to get_local_component is coming from ggf::transfer_from_all at: int const lc2 = h.get_local_component(rl2,c2); where c2 is from psend.component. So it looks like there is an off-by-one error somewhere.

Very strange. This code should be quite solid by now. psend is set in the file dh.cc in thorn Carpet/CarpetLib; there is one (large) routine that calculates the communication schedule. Some of the indexing errors there in the past included confusing the number of components on different refinement levels, which led to indexing errors such as the one you describe.

The bad component numbers are not coming from: preg.component = tmpncomps.AT(m)++; in Carpet/src/Recompose.cc

Where else are the component numbers assigned?

---

Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

---

Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users

Erik,

I've attached my grid_structure output.

It seems to get to: Evolving iteration 3073 at t=12.0039 (as you had suggested, I am using 10 levels).

Thanks again, Hal

On Tue, 2011-09-13 at 18:58 -0400, Erik Schnetter wrote:

After how many iterations does the code abort?

-erik

On Tue, Sep 13, 2011 at 6:45 PM, Erik Schnetter schnetter@cct.lsu.edu wrote:

...

Hal

Were you running with multiple processes?

-erik

On Tue, Sep 13, 2011 at 12:00 PM, Hal Finkel hfinkel@anl.gov wrote:

...

Ian,

I'd not found the problem yet. Hopefully this will be the hint we need to get this fixed.

Thanks again, Hal

On Tue, 2011-09-13 at 15:31 +0100, Ian Hawke wrote:

...

Erik, Hal,

Did you have any luck tracking down this error? I've just come back to this and am seeing the same error message; it appears to arise when two grids on a refined level merge, as in:

[3][0][0]   exterior: [-0.010000,-0.020000,-0.020000] :

[0.000000,0.020000,0.020000] : [0.001250,0.001250,0.001250] [3][0][1] exterior: [0.010000,-0.020000,-0.020000] : [0.040000,0.020000,0.020000] : [0.001250,0.001250,0.001250]

becomes

[3][0][0] exterior: [-0.010000,-0.020000,-0.020000] : [0.040000,0.020000,0.020000] : [0.001250,0.001250,0.001250]

It seems that the old grids are destroyed before the data is copied/populated in Recompose (either that or the old grid structure is not referred to in the data transfer).

Ian

On 07/09/11 01:03, Erik Schnetter wrote:

...

Hal

This is where numbers are assigned to components. The communication schedule decides which component needs to send data to which other component (which may be located on another process or not); this schedule is created for each refinement level independently, and may (if there is an error) refer to component numbers that don't exist. This schedule is set up in dh.cc.

Can you send me the example you are currently running (your source code and parameter file)? I will try to give it a try.

-erik

On Tue, Sep 6, 2011 at 7:44 PM, Hal Finkelhfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 15:16 -0400, Erik Schnetter wrote: > On Thu, Sep 1, 2011 at 3:05 PM, Hal Finkelhfinkel@anl.gov wrote: >> On Thu, 2011-09-01 at 14:25 -0400, Erik Schnetter wrote: >>> On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkelhfinkel@anl.gov wrote: >>>> On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote: >>>>> On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkelhfinkel@anl.gov wrote: >>>>>> On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote: >>>>>>> On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkelhfinkel@anl.gov wrote: >>>>>>>> Could I also decrease the block size? I currently have >>>>>>>> CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? >>>>>>>> Is there a restriction based on the number of ghost points? >>>>>>> Yes, you can reduce the block size. I assume that both the regridding >>>>>>> operation and the time evolution will become slower if you do that, >>>>>>> because more blocks will have to be handled. >>>>>> Regardless of what I do, once we get past the first coarse time step, >>>>>> the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] >>>>>> Regridding map 0...". >>>>>> >>>>>> Overall, it is in dh::regrid(do_init=true). It spends most of its time >>>>>> in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: >>>>>> for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ >>>>>> nsi). The normalize() function does exit, however, so it is not hanging >>>>>> in that function. >>>>>> >>>>>> The core problem seems to be that it takes a long time to execute: >>>>>> boxes = boxes .shift(-dir) - boxes; >>>>>> in dh::regrid(do_init=true). Probably because boxes has 129064 elements. >>>>>> The coarse grid is now only 30^3 and I've left the regrid box size at 4. >>>>>> I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 >>>>>> ~ 420 refinement regions. >>>>>> >>>>>> What is the best way to figure out what is going on? >>>>> Hal >>>>> >>>>> Yes, this function is very slow. I did not expect it to be >>>>> prohibitively slow. Are you compiling with optimisation enabled? >>>> I've tried with optimizations enabled (and without for debugging). >>>> >>>>> The bboxset represents the set of refined regions, and it is >>>>> internally represented as a list of bboxes (regions). Carpet performs >>>>> set operations on these (intersection, union, complement, etc.) to >>>>> determine the communication schedule, i.e. which ghost zones of which >>>>> bbox need to be filled from which other bbox. Unfortunately, the >>>>> algorithm used for this is O(n^2) in the number of refined regions, >>>>> and set operations when implemented via lists themselves are O(n^2) in >>>>> the set size, leading to a rather unfortunate overall complexity. The >>>>> only cure is to reduce the number of bboxes (make them larger) and to >>>>> regrid fewer times. >>>> This is what I suspected, but nevertheless, is there something wrong? >>>> How many boxes do you expect that I should have? The reason that it does >>>> not finish, even with optimizations, is that there are 129K boxes in the >>>> loop (that's at least 16 billion box normalizations?). >>>> >>>> The coarse grid is only 30^3, and the regrid box size is 4, so at >>>> maximum, there should be ~400 level one boxes. Even if some of those >>>> have level 2 boxes, I don't understand how there could be 129K boxes. >>> The refinement structure itself should have one bbox per refined 4^3 >>> box, and both CarpetRegrid2 and CarpetLib would try to combine these >>> into fewer boxes where possible, i.e. where one can form rectangles or >>> larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on >>> level one. >> That makes sense. I think that there is a bug somewhere which is causing >> the box set to be much too big. Furthermore, it does not happen on every >> run, only sometimes. When it does not happen, I hit another bug after a >> few coarse timesteps: >> >> I get a range-check exception from std::vector in a call to: >> gh::get_local_component (rl=1, c=8) >> the problem is that this returns: >> local_components_.AT(rl).AT(c); >> and local_components_[1].size() is 8 >> The call to get_local_component is coming from ggf::transfer_from_all >> at: >> int const lc2 = h.get_local_component(rl2,c2); >> where c2 is from psend.component. >> So it looks like there is an off-by-one error somewhere. > Very strange. This code should be quite solid by now. psend is set in > the file dh.cc in thorn Carpet/CarpetLib; there is one (large) routine > that calculates the communication schedule. Some of the indexing > errors there in the past included confusing the number of components > on different refinement levels, which led to indexing errors such as > the one you describe. The bad component numbers are not coming from: preg.component = tmpncomps.AT(m)++; in Carpet/src/Recompose.cc

Where else are the component numbers assigned?

---

Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

---

Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users

-- Erik Schnetter schnetter@cct.lsu.edu http://www.cct.lsu.edu/~eschnett/

On Wed, Sep 14, 2011 at 9:41 PM, Erik Schnetter schnetter@cct.lsu.edu wrote:

Hal

I believe I found the cause; I am attaching a patch. Could you try the patch and see whether it works? Unfortunately I have other, unrelated changes in my source tree as well.

The problem is the mapping from component number to process number. This mapping may change during regridding. Carpet always used the mapping of the current grid structure, and never used the mapping from the previous grid structure, and this can lead to failure during regridding.

This is related not only to the mapping from components to processes, but also the mapping from global component indices to process-local component indices. While both were (probably) wrong, it is the latter that caused problems for you: The number of components decreased during regridding, and component #8 was looked up in the new mapping instead of the old.

It boggles my mind that such a problem would not manifest itself itself in other simulations as well, either leading to catastrophic errors in simulations or in outright crashes. I assume that the reason is that all other simulations have many fewer components than you, and the mapping is thus trivial and never changes, so that this error is not triggered.

-erik

-- Erik Schnetter schnetter@cct.lsu.edu   http://www.cct.lsu.edu/~eschnett/

On 15/09/11 02:41, Erik Schnetter wrote:

Hal

I believe I found the cause; I am attaching a patch. Could you try the patch and see whether it works? Unfortunately I have other, unrelated changes in my source tree as well.

Fixes the problems in my tests with GRHydro as well (of course, it shows up all the problems in the error measure now...).

Ian

Erik,

The problem occurs with just one process.

Thanks again, Hal

On Tue, 2011-09-13 at 18:45 -0400, Erik Schnetter wrote:

Hal

Were you running with multiple processes?

-erik

On Tue, Sep 13, 2011 at 12:00 PM, Hal Finkel hfinkel@anl.gov wrote:

...

Ian,

I'd not found the problem yet. Hopefully this will be the hint we need to get this fixed.

Thanks again, Hal

On Tue, 2011-09-13 at 15:31 +0100, Ian Hawke wrote:

...

Erik, Hal,

Did you have any luck tracking down this error? I've just come back to this and am seeing the same error message; it appears to arise when two grids on a refined level merge, as in:

[3][0][0]   exterior: [-0.010000,-0.020000,-0.020000] :

[0.000000,0.020000,0.020000] : [0.001250,0.001250,0.001250] [3][0][1] exterior: [0.010000,-0.020000,-0.020000] : [0.040000,0.020000,0.020000] : [0.001250,0.001250,0.001250]

becomes

[3][0][0] exterior: [-0.010000,-0.020000,-0.020000] : [0.040000,0.020000,0.020000] : [0.001250,0.001250,0.001250]

It seems that the old grids are destroyed before the data is copied/populated in Recompose (either that or the old grid structure is not referred to in the data transfer).

Ian

On 07/09/11 01:03, Erik Schnetter wrote:

...

Hal

This is where numbers are assigned to components. The communication schedule decides which component needs to send data to which other component (which may be located on another process or not); this schedule is created for each refinement level independently, and may (if there is an error) refer to component numbers that don't exist. This schedule is set up in dh.cc.

Can you send me the example you are currently running (your source code and parameter file)? I will try to give it a try.

-erik

On Tue, Sep 6, 2011 at 7:44 PM, Hal Finkelhfinkel@anl.gov wrote:

...

On Thu, 2011-09-01 at 15:16 -0400, Erik Schnetter wrote:

...

On Thu, Sep 1, 2011 at 3:05 PM, Hal Finkelhfinkel@anl.gov wrote: > On Thu, 2011-09-01 at 14:25 -0400, Erik Schnetter wrote: >> On Thu, Sep 1, 2011 at 11:51 AM, Hal Finkelhfinkel@anl.gov wrote: >>> On Thu, 2011-09-01 at 11:37 -0400, Erik Schnetter wrote: >>>> On Thu, Sep 1, 2011 at 10:53 AM, Hal Finkelhfinkel@anl.gov wrote: >>>>> On Tue, 2011-08-30 at 21:06 -0400, Erik Schnetter wrote: >>>>>> On Tue, Aug 30, 2011 at 5:28 PM, Hal Finkelhfinkel@anl.gov wrote: >>>>>>> Could I also decrease the block size? I currently have >>>>>>> CarpetRegrid2::adaptive_block_size = 4, could it be smaller than that? >>>>>>> Is there a restriction based on the number of ghost points? >>>>>> Yes, you can reduce the block size. I assume that both the regridding >>>>>> operation and the time evolution will become slower if you do that, >>>>>> because more blocks will have to be handled. >>>>> Regardless of what I do, once we get past the first coarse time step, >>>>> the program seems to "hang" at "INFO (Carpet): [ml=0][rl=0][m=0][tl=0] >>>>> Regridding map 0...". >>>>> >>>>> Overall, it is in dh::regrid(do_init=true). It spends most of its time >>>>> in bboxset<int, 3>::normalize() and, specifically, mostly in the loop: >>>>> for (typename bset::iterator nsi = nbs.begin(); nsi != nbs.end(); ++ >>>>> nsi). The normalize() function does exit, however, so it is not hanging >>>>> in that function. >>>>> >>>>> The core problem seems to be that it takes a long time to execute: >>>>> boxes = boxes .shift(-dir) - boxes; >>>>> in dh::regrid(do_init=true). Probably because boxes has 129064 elements. >>>>> The coarse grid is now only 30^3 and I've left the regrid box size at 4. >>>>> I'd think, then, that the coarse grid should have a maximum of 30^3/4^3 >>>>> ~ 420 refinement regions. >>>>> >>>>> What is the best way to figure out what is going on? >>>> Hal >>>> >>>> Yes, this function is very slow. I did not expect it to be >>>> prohibitively slow. Are you compiling with optimisation enabled? >>> I've tried with optimizations enabled (and without for debugging). >>> >>>> The bboxset represents the set of refined regions, and it is >>>> internally represented as a list of bboxes (regions). Carpet performs >>>> set operations on these (intersection, union, complement, etc.) to >>>> determine the communication schedule, i.e. which ghost zones of which >>>> bbox need to be filled from which other bbox. Unfortunately, the >>>> algorithm used for this is O(n^2) in the number of refined regions, >>>> and set operations when implemented via lists themselves are O(n^2) in >>>> the set size, leading to a rather unfortunate overall complexity. The >>>> only cure is to reduce the number of bboxes (make them larger) and to >>>> regrid fewer times. >>> This is what I suspected, but nevertheless, is there something wrong? >>> How many boxes do you expect that I should have? The reason that it does >>> not finish, even with optimizations, is that there are 129K boxes in the >>> loop (that's at least 16 billion box normalizations?). >>> >>> The coarse grid is only 30^3, and the regrid box size is 4, so at >>> maximum, there should be ~400 level one boxes. Even if some of those >>> have level 2 boxes, I don't understand how there could be 129K boxes. >> The refinement structure itself should have one bbox per refined 4^3 >> box, and both CarpetRegrid2 and CarpetLib would try to combine these >> into fewer boxes where possible, i.e. where one can form rectangles or >> larger cubes. I would thus expect no more than (30/4)^2 = 64 bboxes on >> level one. > That makes sense. I think that there is a bug somewhere which is causing > the box set to be much too big. Furthermore, it does not happen on every > run, only sometimes. When it does not happen, I hit another bug after a > few coarse timesteps: > > I get a range-check exception from std::vector in a call to: > gh::get_local_component (rl=1, c=8) > the problem is that this returns: > local_components_.AT(rl).AT(c); > and local_components_[1].size() is 8 > The call to get_local_component is coming from ggf::transfer_from_all > at: > int const lc2 = h.get_local_component(rl2,c2); > where c2 is from psend.component. > So it looks like there is an off-by-one error somewhere. Very strange. This code should be quite solid by now. psend is set in the file dh.cc in thorn Carpet/CarpetLib; there is one (large) routine that calculates the communication schedule. Some of the indexing errors there in the past included confusing the number of components on different refinement levels, which led to indexing errors such as the one you describe.

The bad component numbers are not coming from: preg.component = tmpncomps.AT(m)++; in Carpet/src/Recompose.cc

Where else are the component numbers assigned?

---

Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 1-630-252-0023 hfinkel@anl.gov

---

Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users

Hal

This may mean that level_mask does not have storage at all times. CarpetRegrid2 does not allocate storage for this variable -- maybe it should, when CarpetRegrid2::adaptive_refinement is set to yes?

-erik

On Wed, Dec 14, 2011 at 4:07 PM, Hal Finkel hfinkel@anl.gov wrote:

On Thu, 2011-12-08 at 22:38 +0100, Eloisa Bentivegna wrote:

...

Hello,

I'd like to revive a discussion that went over the list a few months

ago, regarding Carpet's AMR capability.

...

I understand that the refinement has to be triggered through the

level_mask grid function. The example that Hal gave is below.

...

On Aug 30, 2011, at 9:37 PM, Hal Finkel wrote:

...

lvlmskSet = { Name -> "SFBubble_SetLevelMask", Schedule -> {"AT preregrid"}, Where -> Interior, Shorthands -> { adchix, adchiy, adchiz }, Equations -> { adchix -> fabs[dx PD[chi,1]], adchiy -> fabs[dy PD[chi,2]], adchiz -> fabs[dz PD[chi,3]],

           "level_mask" -> (
                   Max[
                           adchix, adchiy, adchiz
                   ]/dchimax
           ) /. Max[a_, b_, c_] -> fmax[a, fmax[b, c]]
   }

};

FWIW, it looks like some time over the last few months this stopped working: setting level_mask in this way does not lead to regridding (even if it is > n in some places), and outputting the result using CarpetIOScalar has an odd result:

carpetregrid2::level_mask.maximum.asc starts with: 0 0 -1.79769313486232e+308 1 1 1.30685299729774 2 2 1.31094960089828 ...

carpetregrid2::level_mask.minimum.asc starts with: 0 0 1.79769313486232e+308 1 1 0 2 2 0 ...

carpetregrid2::level_mask.average.asc has: 0 0 -nan 1 1 -nan 2 2 -nan ...

What is going on here?

It might also be worth noting, that to get level_mask to work with Kranc in the above example, I added the following after the call to CreateKrancThornTT:

Module[{fp = OpenAppend[dirname <> "/schedule.ccl"]}, WriteString[fp, "[NewLine]" <> "STORAGE: CarpetRegrid2::level_mask"]; Close[fp] ];

-Hal

...

My questions are:

1. Is preregrid the right location to set this function?
2. Does this suffice? I assume that once regridding is triggered, the

new levels will have to be populated (via interpolation?). Is level_mask treated like all other grid functions?

...

3. In the case above, the loop only goes over interior points, since the

mask is set to be the derivative of a grid function. This implies that a sync is necessary before the mask can be used for regridding. Is the user thorn a good place to request this?

...

Thanks, Eloisa

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory

---

Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users

It should suffice if this is a global storage allocation, i.e. if it is outside of any schedule statement.

-erik

On Wed, Dec 14, 2011 at 4:58 PM, Hal Finkel hfinkel@anl.gov wrote:

On Wed, 2011-12-14 at 16:33 -0500, Erik Schnetter wrote:

...

Hal

This may mean that level_mask does not have storage at all times. CarpetRegrid2 does not allocate storage for this variable -- maybe it should, when CarpetRegrid2::adaptive_refinement is set to yes?

I agree that it should do that; but I am allocating storage for the variable right now (as I noted at the end of the e-mail) by adding to the end of my schedule file: STORAGE: CarpetRegrid2::level_mask

Is that not enough?

Thanks again, Hal

...

-erik

On Wed, Dec 14, 2011 at 4:07 PM, Hal Finkel hfinkel@anl.gov wrote: On Thu, 2011-12-08 at 22:38 +0100, Eloisa Bentivegna wrote: > Hello, > > I'd like to revive a discussion that went over the list a few months ago, regarding Carpet's AMR capability. > > I understand that the refinement has to be triggered through the level_mask grid function. The example that Hal gave is below. > > On Aug 30, 2011, at 9:37 PM, Hal Finkel wrote: > > > lvlmskSet = > > { > > Name -> "SFBubble_SetLevelMask", > > Schedule -> {"AT preregrid"}, > > Where -> Interior, > > Shorthands -> { adchix, adchiy, adchiz }, > > Equations -> > > { > > adchix -> fabs[dx PD[chi,1]], > > adchiy -> fabs[dy PD[chi,2]], > > adchiz -> fabs[dz PD[chi,3]], > > > > "level_mask" -> ( > > Max[ > > adchix, adchiy, adchiz > > ]/dchimax > > ) /. Max[a_, b_, c_] -> fmax[a, fmax[b, c]] > > } > > };

    FWIW, it looks like some time over the last few months this
    stopped
    working: setting level_mask in this way does not lead to
    regridding
    (even if it is > n in some places), and outputting the result
    using
    CarpetIOScalar has an odd result:

    carpetregrid2::level_mask.maximum.asc starts with:
    0 0 -1.79769313486232e+308
    1 1 1.30685299729774
    2 2 1.31094960089828
    ...

    carpetregrid2::level_mask.minimum.asc starts with:
    0 0 1.79769313486232e+308
    1 1 0
    2 2 0
    ...

    carpetregrid2::level_mask.average.asc has:
    0 0 -nan
    1 1 -nan
    2 2 -nan
    ...

    What is going on here?

    It might also be worth noting, that to get level_mask to work
    with Kranc
    in the above example, I added the following after the call to
    CreateKrancThornTT:

    Module[{fp = OpenAppend[dirname <> "/schedule.ccl"]},
           WriteString[fp, "\[NewLine]" <> "STORAGE:
    CarpetRegrid2::level_mask"];
           Close[fp]
    ];

     -Hal

    >
    > My questions are:
    >
    > 1) Is preregrid the right location to set this function?
    > 2) Does this suffice? I assume that once regridding is
    triggered, the new levels will have to be populated (via
    interpolation?). Is level_mask treated like all other grid
    functions?
    > 3) In the case above, the loop only goes over interior
    points, since the mask is set to be the derivative of a grid
    function. This implies that a sync is necessary before the
    mask can be used for regridding. Is the user thorn a good
    place to request this?
    >
    > Thanks,
    > Eloisa


    --
    Hal Finkel
    Postdoctoral Appointee
    Leadership Computing Facility
    Argonne National Laboratory


    _______________________________________________
    Users mailing list
    Users@einsteintoolkit.org
    http://lists.einsteintoolkit.org/mailman/listinfo/users

-- Erik Schnetter schnetter@cct.lsu.edu http://www.cct.lsu.edu/~eschnett/

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory

On Wed, 2011-12-14 at 18:27 -0500, Erik Schnetter wrote:

You are not setting any values on the boundary. Is that intentional?

Currently, this is because I am using neighboring values in the calculation, so I can't do that on the boundary. Should I do that some other way?

Can you look at the mask everywhere, to ensure that no value is a nan?

There are NaNs on the grid. They exist whereever any coordinate has its highest value (not for 0, just for the highest value). Looks almost like an off-by-one error somewhere.

The max and min intrinsics may not pick up nans.

It seems that they don't ;)

To get things working for initial data, you may need to set the mask in the preregridinitial bin as well, or you need to ensure that carpet doesn't use the AMR mechanism at this time, or that another refinement mechanism sets up a sufficiently refined grid structure. (Once the grid structure is too coarse, things go wrong, and calculating derivative may lead to nans, so that level_mask cannot be calculated any more...)

Currently, no regridding is happening at all, but this may be a function of the NaNs.

Maybe you need a secondary mechanism while setting level_mask, to ensure there is a minimum amount of refinement present around certain features? I can't tell, this depends on your physical system.

The NaNs are not, as far as I can tell, coming from a calculation, but rather they're coming from uninitialized data.

Thanks again, Hal

-erik

On Wed, Dec 14, 2011 at 4:07 PM, Hal Finkel hfinkel@anl.gov wrote: On Thu, 2011-12-08 at 22:38 +0100, Eloisa Bentivegna wrote: > Hello, > > I'd like to revive a discussion that went over the list a few months ago, regarding Carpet's AMR capability. > > I understand that the refinement has to be triggered through the level_mask grid function. The example that Hal gave is below. > > On Aug 30, 2011, at 9:37 PM, Hal Finkel wrote: > > > lvlmskSet = > > { > > Name -> "SFBubble_SetLevelMask", > > Schedule -> {"AT preregrid"}, > > Where -> Interior, > > Shorthands -> { adchix, adchiy, adchiz }, > > Equations -> > > { > > adchix -> fabs[dx PD[chi,1]], > > adchiy -> fabs[dy PD[chi,2]], > > adchiz -> fabs[dz PD[chi,3]], > > > > "level_mask" -> ( > > Max[ > > adchix, adchiy, adchiz > > ]/dchimax > > ) /. Max[a_, b_, c_] -> fmax[a, fmax[b, c]] > > } > > };

    FWIW, it looks like some time over the last few months this
    stopped
    working: setting level_mask in this way does not lead to
    regridding
    (even if it is > n in some places), and outputting the result
    using
    CarpetIOScalar has an odd result:

    carpetregrid2::level_mask.maximum.asc starts with:
    0 0 -1.79769313486232e+308
    1 1 1.30685299729774
    2 2 1.31094960089828
    ...

    carpetregrid2::level_mask.minimum.asc starts with:
    0 0 1.79769313486232e+308
    1 1 0
    2 2 0
    ...

    carpetregrid2::level_mask.average.asc has:
    0 0 -nan
    1 1 -nan
    2 2 -nan
    ...

    What is going on here?

    It might also be worth noting, that to get level_mask to work
    with Kranc
    in the above example, I added the following after the call to
    CreateKrancThornTT:

    Module[{fp = OpenAppend[dirname <> "/schedule.ccl"]},
           WriteString[fp, "\[NewLine]" <> "STORAGE:
    CarpetRegrid2::level_mask"];
           Close[fp]
    ];

     -Hal

    >
    > My questions are:
    >
    > 1) Is preregrid the right location to set this function?
    > 2) Does this suffice? I assume that once regridding is
    triggered, the new levels will have to be populated (via
    interpolation?). Is level_mask treated like all other grid
    functions?
    > 3) In the case above, the loop only goes over interior
    points, since the mask is set to be the derivative of a grid
    function. This implies that a sync is necessary before the
    mask can be used for regridding. Is the user thorn a good
    place to request this?
    >
    > Thanks,
    > Eloisa


    --
    Hal Finkel
    Postdoctoral Appointee
    Leadership Computing Facility
    Argonne National Laboratory


    _______________________________________________
    Users mailing list
    Users@einsteintoolkit.org
    http://lists.einsteintoolkit.org/mailman/listinfo/users

-- Erik Schnetter schnetter@cct.lsu.edu http://www.cct.lsu.edu/~eschnett/

On Wed, 2011-12-14 at 19:48 -0500, Erik Schnetter wrote:

On Wed, Dec 14, 2011 at 6:53 PM, Hal Finkel hfinkel@anl.gov wrote: On Wed, 2011-12-14 at 18:27 -0500, Erik Schnetter wrote: > You are not setting any values on the boundary. Is that intentional?

    Currently, this is because I am using neighboring values in
    the
    calculation, so I can't do that on the boundary. Should I do
    that some
    other way?

Setting the boundary to zero should be good enough. (The boundary values should not be used -- but I don't recall whether this is the case.)

Maybe this is another problem with the periodic boundary conditions? The problem seems to appear in other fields too. I've attached some images from my test problem. One shows the field configuration (this one looks like a ring -- it is a slice through a bubble). The second one shows the T_00 computed from that. As you can see, the values near the extremal indicies are wrong. At the next (half) time step, looking at the field data from level 1, the same phenomonon can be seen in the level 1 boxes.

These were all done in Kranc, so I did not do any real coding myself ;)

What do you think?

Thanks again, Hal

-erik

-- Erik Schnetter schnetter@cct.lsu.edu http://www.cct.lsu.edu/~eschnett/

Eloisa

Thanks for digging into this.

I didn't think of syncing, but yes, syncing would be necessary, and so would be setting all boundaries. In your case, this would be applying periodic boundary conditions only (and not setting anything to zero since you don't have an outer boundary). Of course, you need to use sufficiently many processes for this because of the bug in thorn Periodic.

-erik

On Thu, Dec 15, 2011 at 9:39 AM, Eloisa Bentivegna bentivegna@cct.lsu.eduwrote:

On Dec 15, 2011, at 5:58 AM, Hal Finkel wrote:

...

On Wed, 2011-12-14 at 19:48 -0500, Erik Schnetter wrote:

...

On Wed, Dec 14, 2011 at 6:53 PM, Hal Finkel hfinkel@anl.gov wrote: On Wed, 2011-12-14 at 18:27 -0500, Erik Schnetter wrote:

...

You are not setting any values on the boundary. Is that

   intentional?


   Currently, this is because I am using neighboring values in
   the
   calculation, so I can't do that on the boundary. Should I do
   that some
   other way?

Setting the boundary to zero should be good enough. (The boundary values should not be used -- but I don't recall whether this is the case.)

Maybe this is another problem with the periodic boundary conditions? The problem seems to appear in other fields too. I've attached some images from my test problem. One shows the field configuration (this one looks like a ring -- it is a slice through a bubble). The second one shows the T_00 computed from that. As you can see, the values near the extremal indicies are wrong. At the next (half) time step, looking at the field data from level 1, the same phenomonon can be seen in the level 1 boxes.

These were all done in Kranc, so I did not do any real coding myself ;)

What do you think?

Hi all!

I believe that what Hal is doing with the allocation of level_mask is correct; what seems to be problematic is setting its value.

The problem is that this function cannot be calculated pointwise, and the way it's currently set leaves some parts of the grid uninitialized, which then leads to poison and ultimately to all sorts of weird behavior including nans at the boundaries, lack of regridding when expected, and so on. I wouldn't pay too much attention to all these symptoms (although it would be nice of course for the AMR logic to detect that the mask is nan and issue an error message). The origin is in the incomplete initialization of the mask.

To convince myself of this, I switched poisoning off, and observed no odd behavior. I then switched it back on and ran a number of grid configurations where the mask was set pointwise (say, to coincide with one of the evolution variables), and again I had no trouble. As for the actual mask (based on the derivative of a grid function), I played a long time with scheduling the filling of interior points and boundaries, but always observed the issue that Hal is reporting. Ideally, I'd think that filling the interior, syncing, and finally applying outer/symmetry boundary conditions would work, but that doesn't seem to be the case.

Erik: when you suggest to set the boundary explicitly, do you mean the outer boundary? Since in this case we only have a symmetry boundary (periodic), do you mean we should populate that part of the grid independently of the symmetry thorn?

Thanks, Eloisa

Hal

You should be able to apply boundary conditions (i.e. synchronised, apply symmetry conditions) in your code. You will also need to use sufficiently many processes to ensure that periodicity is applied correctly; this circumvents the bug/missing feature in the Periodic boundary conditions.

-erik

On Thu, Dec 15, 2011 at 3:14 PM, Hal Finkel hfinkel@anl.gov wrote:

On Thu, 2011-12-15 at 11:17 -0500, Erik Schnetter wrote:

...

Eloisa

Thanks for digging into this.

I didn't think of syncing, but yes, syncing would be necessary, and so would be setting all boundaries. In your case, this would be applying periodic boundary conditions only (and not setting anything to zero since you don't have an outer boundary). Of course, you need to use sufficiently many processes for this because of the bug in thorn Periodic.

Erik, Eloisa, et al.,

Thank you very much for looking at this. Just so that I'm clear, is there a way that I can get this to work now, or am I waiting on another bug fix?

-Hal

...

-erik

On Thu, Dec 15, 2011 at 9:39 AM, Eloisa Bentivegna bentivegna@cct.lsu.edu wrote: On Dec 15, 2011, at 5:58 AM, Hal Finkel wrote:

    > On Wed, 2011-12-14 at 19:48 -0500, Erik Schnetter wrote:
    >> On Wed, Dec 14, 2011 at 6:53 PM, Hal Finkel
    <hfinkel@anl.gov> wrote:
    >>        On Wed, 2011-12-14 at 18:27 -0500, Erik Schnetter
    wrote:
    >>> You are not setting any values on the boundary. Is that
    >>        intentional?
    >>
    >>
    >>        Currently, this is because I am using neighboring
    values in
    >>        the
    >>        calculation, so I can't do that on the boundary.
    Should I do
    >>        that some
    >>        other way?
    >>
    >>
    >> Setting the boundary to zero should be good enough. (The
    boundary
    >> values should not be used -- but I don't recall whether
    this is the
    >> case.)
    >
    > Maybe this is another problem with the periodic boundary
    conditions?
    > The problem seems to appear in other fields too. I've
    attached some
    > images from my test problem. One shows the field
    configuration (this one
    > looks like a ring -- it is a slice through  a bubble). The
    second one
    > shows the T_00 computed from that. As you can see, the
    values near the
    > extremal indicies are wrong. At the next (half) time step,
    looking at
    > the field data from level 1, the same phenomonon can be seen
    in the
    > level 1 boxes.
    >
    > These were all done in Kranc, so I did not do any real
    coding myself ;)
    >
    > What do you think?


    Hi all!

    I believe that what Hal is doing with the allocation of
    level_mask is correct; what seems to be problematic is setting
    its value.

    The problem is that this function cannot be calculated
    pointwise, and the way it's currently set leaves some parts of
    the grid uninitialized, which then leads to poison and
    ultimately to all sorts of weird behavior including nans at
    the boundaries, lack of regridding when expected, and so on. I
    wouldn't pay too much attention to all these symptoms
    (although it would be nice of course for the AMR logic to
    detect that the mask is nan and issue an error message). The
    origin is in the incomplete initialization of the mask.

    To convince myself of this, I switched poisoning off, and
    observed no odd behavior. I then switched it back on and ran a
    number of grid configurations where the mask was set pointwise
    (say, to coincide with one of the evolution variables), and
    again I had no trouble. As for the actual mask (based on the
    derivative of a grid function), I played a long time with
    scheduling the filling of interior points and boundaries, but
    always observed the issue that Hal is reporting. Ideally, I'd
    think that filling the interior, syncing, and finally applying
    outer/symmetry boundary conditions would work, but that
    doesn't seem to be the case.

    Erik: when you suggest to set the boundary explicitly, do you
    mean the outer boundary? Since in this case we only have a
    symmetry boundary (periodic), do you mean we should populate
    that part of the grid independently of the symmetry thorn?

    Thanks,
    Eloisa

-- Erik Schnetter schnetter@cct.lsu.edu http://www.cct.lsu.edu/~eschnett/

-- Hal Finkel Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory

On Dec 15, 2011, at 9:50 PM, Hal Finkel wrote:

Can you be more specific?

I think that what Erik is suggesting is to fill all the grid parts that are not covered by the SetLevelMask calculation. These entail:

1) The interprocessor and mesh refinement boundaries, which are filled by requiring that level_mask is synchronised upon exit from SetLevelMask (basically this is done by adding "SYNC: level_mask" to the schedule block for SetLevelMask in schedule.ccl); 2) The symmetry boundaries, which are filled by applying the corresponding conditions right after SetLevelMask (this is accomplished by scheduling the SelectBoundConds function and ApplyBCs group, in this order, after SetLevelMask).

It's worth stressing that this approach doesn't currently work for me; I plan on taking a closer look in the next few days.

And is a sufficient number just to make sure that no one processor holds to opposing faces of the cube?

Yes. I'm not sure whether this number can be calculated a priori, especially in the case where the grid is worked out at runtime and may be quite dynamical. However, I'm reviewing Erik's patch introducing some login in Periodic to check exactly this. With this patch, you'll know (a posteriori) whether you have a sufficient number of processors.

Thanks, Eloisa