Hi,
What is the current situation with evolving grid arrays with MoL? Over the past 7 years, this topic has come up repeatedly on the mailing lists, and some of the old information is probably out of date. Is there a way to do this without using any hacks, or are hacks still necessary?
Is there an example somewhere of how to do this, or does someone have an example they would be willing to share?
Hi Ian,
As far as I know this wors fine. I have used it in EHFinder to evolve the generators of the horizon, where the coordinates are stored as 1-d grid arrays (like xg and its rhs dxg). However, I haven't run that code in a while, so I can't guarantee that it still works. I'm also using it in a thorn to integrate geodesics in Schwarzschild (for the self-force project). This thorn does work currently but is not yet in the EinsteinToolkit. If you want to take a look at it, just let me know.
Cheers,
Peter
On Mon, 23 Sep 2013, Ian Hinder wrote:
Hi, What is the current situation with evolving grid arrays with MoL? Over the past 7 years, this topic has come up repeatedly on the mailing lists, and some of the old information is probably out of date. Is there a way to do this without using any hacks, or are hacks still necessary?
Is there an example somewhere of how to do this, or does someone have an example they would be willing to share?
-- Ian Hinder http://numrel.aei.mpg.de/people/hinder
does it work with AMR and MOL? We use GA's to evolve the puncture locations. There are awful hacks there to make sure the GAs are evolved correctly with MOL.
On 09/23/2013 10:51 AM, Peter Diener wrote:
Hi Ian,
As far as I know this wors fine. I have used it in EHFinder to evolve the generators of the horizon, where the coordinates are stored as 1-d grid arrays (like xg and its rhs dxg). However, I haven't run that code in a while, so I can't guarantee that it still works. I'm also using it in a thorn to integrate geodesics in Schwarzschild (for the self-force project). This thorn does work currently but is not yet in the EinsteinToolkit. If you want to take a look at it, just let me know.
Cheers,
Peter
On Mon, 23 Sep 2013, Ian Hinder wrote:
Hi, What is the current situation with evolving grid arrays with MoL? Over the past 7 years, this topic has come up repeatedly on the mailing lists, and some of the old information is probably out of date. Is there a way to do this without using any hacks, or are hacks still necessary?
Is there an example somewhere of how to do this, or does someone have an example they would be willing to share?
-- Ian Hinder http://numrel.aei.mpg.de/people/hinder
Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users
Hi,
No, you're right. The examples I mentioned are using either uni-grid or multi-patch without AMR. I didn't think of the complication of adding AMR to the mix.
Cheers,
Peter
On Mon, 23 Sep 2013, Yosef Zlochower wrote:
does it work with AMR and MOL? We use GA's to evolve the puncture locations. There are awful hacks there to make sure the GAs are evolved correctly with MOL.
On 09/23/2013 10:51 AM, Peter Diener wrote:
Hi Ian,
As far as I know this wors fine. I have used it in EHFinder to evolve the generators of the horizon, where the coordinates are stored as 1-d grid arrays (like xg and its rhs dxg). However, I haven't run that code in a while, so I can't guarantee that it still works. I'm also using it in a thorn to integrate geodesics in Schwarzschild (for the self-force project). This thorn does work currently but is not yet in the EinsteinToolkit. If you want to take a look at it, just let me know.
Cheers,
Peter
On Mon, 23 Sep 2013, Ian Hinder wrote:
Hi, What is the current situation with evolving grid arrays with MoL? Over the past 7 years, this topic has come up repeatedly on the mailing lists, and some of the old information is probably out of date. Is there a way to do this without using any hacks, or are hacks still necessary?
Is there an example somewhere of how to do this, or does someone have an example they would be willing to share?
-- Ian Hinder http://numrel.aei.mpg.de/people/hinder
Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users
-- Dr. Yosef Zlochower Center for Computational Relativity and Gravitation Associate Professor School of Mathematical Sciences Rochester Institute of Technology 85 Lomb Memorial Drive Rochester, NY 14623
Office:74-2067 Phone: +1 585-475-6103
yosef@astro.rit.edu
CONFIDENTIALITY NOTE: The information transmitted, including attachments, is intended only for the person(s) or entity to which it is addressed and may contain confidential and/or privileged material. Any review, retransmission, dissemination or other use of, or taking of any action in reliance upon this information by persons or entities other than the intended recipient is prohibited. If you received this in error, please contact the sender and destroy any copies of this information.
On 23 Sep 2013, at 18:06, Peter Diener diener@cct.lsu.edu wrote:
Hi,
No, you're right. The examples I mentioned are using either uni-grid or multi-patch without AMR. I didn't think of the complication of adding AMR to the mix.
So unless anything has changed, the latest information I have is from Erik's email from September 2007:
On 28 Sep 2007, at 18:10, Erik Schnetter schnetter@cct.lsu.edu wrote:
I am now evolving geodesics with the following work-around:
I use MoL. Each time the RHS is calculated, I check whether the current refinement level is the finest level. If so, the RHS is calculated, otherwise the RHS is set to zero.
The current refinement level is accessed by
interface.ccl: USES INCLUDE: carpet.hh
#include <carpet.hh>
if (Carpet::reflevel == Carpet::reflevels - 1) // do calculation
This works only from C++. The RHS calculating routine has to be scheduled in level mode.
However, I get stuck well before this stage. I try to register the array variable with MoL, but it complains that it is not an array.
Yes, this (i.e. the work-around of scheduling the RHS in level mode, and then only actually executing the routine if this corresponds to global mode) is how I made things work. I haven't tried this recently, though. Maybe a test case would be in order...
Rather straightforward changes to MoL should allow scheduling the RHS in global mode as well, which would then work elegantly. This would probably require to duplicate MoL's integration logic, so that grid functions are integrated in level mode while grid arrays are integrated in global mode.
-erik
On 2013-09-24, at 5:34 , Ian Hinder ian.hinder@aei.mpg.de wrote:
On 23 Sep 2013, at 18:06, Peter Diener diener@cct.lsu.edu wrote:
Hi,
No, you're right. The examples I mentioned are using either uni-grid or multi-patch without AMR. I didn't think of the complication of adding AMR to the mix.
So unless anything has changed, the latest information I have is from Erik's email from September 2007:
On 28 Sep 2007, at 18:10, Erik Schnetter schnetter@cct.lsu.edu wrote:
I am now evolving geodesics with the following work-around:
I use MoL. Each time the RHS is calculated, I check whether the current refinement level is the finest level. If so, the RHS is calculated, otherwise the RHS is set to zero.
The current refinement level is accessed by
interface.ccl: USES INCLUDE: carpet.hh
#include <carpet.hh>
if (Carpet::reflevel == Carpet::reflevels - 1) // do calculation
This works only from C++. The RHS calculating routine has to be scheduled in level mode.
However, I get stuck well before this stage. I try to register the array variable with MoL, but it complains that it is not an array.
-- Ian Hinder http://numrel.aei.mpg.de/people/hinder
Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users
On 28 Sep 2013, at 16:17, Erik Schnetter schnetter@gmail.com wrote:
Yes, this (i.e. the work-around of scheduling the RHS in level mode, and then only actually executing the routine if this corresponds to global mode) is how I made things work. I haven't tried this recently, though. Maybe a test case would be in order...
Rather straightforward changes to MoL should allow scheduling the RHS in global mode as well, which would then work elegantly. This would probably require to duplicate MoL's integration logic, so that grid functions are integrated in level mode while grid arrays are integrated in global mode.
I suspect that the error I ran into during registration was caused by the fact that I was trying to register a vector grid array (or component). If I use a normal (non-vector) grid array, maybe I will get past this first hurdle. Ultimately though, MoL shouldn't care what type of array it is working on.
-erik
On 2013-09-24, at 5:34 , Ian Hinder ian.hinder@aei.mpg.de wrote:
On 23 Sep 2013, at 18:06, Peter Diener diener@cct.lsu.edu wrote:
Hi,
No, you're right. The examples I mentioned are using either uni-grid or multi-patch without AMR. I didn't think of the complication of adding AMR to the mix.
So unless anything has changed, the latest information I have is from Erik's email from September 2007:
On 28 Sep 2007, at 18:10, Erik Schnetter schnetter@cct.lsu.edu wrote:
I am now evolving geodesics with the following work-around:
I use MoL. Each time the RHS is calculated, I check whether the current refinement level is the finest level. If so, the RHS is calculated, otherwise the RHS is set to zero.
The current refinement level is accessed by
interface.ccl: USES INCLUDE: carpet.hh
#include <carpet.hh>
if (Carpet::reflevel == Carpet::reflevels - 1) // do calculation
This works only from C++. The RHS calculating routine has to be scheduled in level mode.
However, I get stuck well before this stage. I try to register the array variable with MoL, but it complains that it is not an array.
-- Ian Hinder http://numrel.aei.mpg.de/people/hinder
Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users
-- Erik Schnetter schnetter@gmail.com 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/.
On 28 Sep 2013, at 16:17, Erik Schnetter schnetter@gmail.com wrote:
Yes, this (i.e. the work-around of scheduling the RHS in level mode, and then only actually executing the routine if this corresponds to global mode) is how I made things work. I haven't tried this recently, though. Maybe a test case would be in order...
Rather straightforward changes to MoL should allow scheduling the RHS in global mode as well, which would then work elegantly. This would probably require to duplicate MoL's integration logic, so that grid functions are integrated in level mode while grid arrays are integrated in global mode.
Is the problem that the grid arrays are not available in global mode, but they are available in level mode? Or is there some other reason that global mode cannot be used? When I try to call the interpolator in MoL_CalcRHS in level mode, I get
Assertion failed: (not (not want_global_mode and num_time_derivs.AT(m) == 0 and need_time_interp.AT(m))), function interpolate_components, file /Users/ian/Cactus/arrangements/Carpet/CarpetInterp/src/interp.cc, line 1353.
This is with InterpNumTimelevels = 1 in the interpolator parameters. I dug around in interp.cc, and came across the interpolation parameter want_global_mode. Setting this to 1 seems to allow the interpolation to succeed, but I don't know if this will break in certain situations. What does this interpolator parameter do?
My problem with registration turned out to be that I was trying to register MoL for grid SCALARs, not grid ARRAYs. A vector of grid arrays seems to work.
On 2013-10-07, at 12:47 , Ian Hinder ian.hinder@aei.mpg.de wrote:
On 28 Sep 2013, at 16:17, Erik Schnetter schnetter@gmail.com wrote:
Yes, this (i.e. the work-around of scheduling the RHS in level mode, and then only actually executing the routine if this corresponds to global mode) is how I made things work. I haven't tried this recently, though. Maybe a test case would be in order...
Rather straightforward changes to MoL should allow scheduling the RHS in global mode as well, which would then work elegantly. This would probably require to duplicate MoL's integration logic, so that grid functions are integrated in level mode while grid arrays are integrated in global mode.
Is the problem that the grid arrays are not available in global mode, but they are available in level mode? Or is there some other reason that global mode cannot be used? When I try to call the interpolator in MoL_CalcRHS in level mode, I get
Assertion failed: (not (not want_global_mode and num_time_derivs.AT(m) == 0 and need_time_interp.AT(m))), function interpolate_components, file /Users/ian/Cactus/arrangements/Carpet/CarpetInterp/src/interp.cc, line 1353.
This is with InterpNumTimelevels = 1 in the interpolator parameters. I dug around in interp.cc, and came across the interpolation parameter want_global_mode. Setting this to 1 seems to allow the interpolation to succeed, but I don't know if this will break in certain situations. What does this interpolator parameter do?
My problem with registration turned out to be that I was trying to register MoL for grid SCALARs, not grid ARRAYs. A vector of grid arrays seems to work.
Grid arrays are available in global mode (and also in level, singlemap, and local mode).
Global mode cannot be used because MoL will still try to evolve them in level mode, i.e. more often. One needs to ensure that grid arrays have a non-zero RHS exactly on the finest grid. Global mode may be equivalent to the coarsest grid, which would then evolve with the wrong step size.
Yes, one cannot interpolate in level mode if this requires data from coarser levels. This is just what level mode means -- access data from the current refinement level only.
Is there a mode "global-finest"? This should work. Or maybe "global-late" for evolution, since finest grid are evolved last? Yes, I think this should work.
Grid scalars are just a funny name for 0-dimensional grid arrays. If MoL doesn't support this, then it would be trivial to add -- the only difference is that they are called "scalar" instead of "array".
-erik
On Mon, Oct 07, 2013 at 01:34:15PM -0400, Erik Schnetter wrote:
Is there a mode "global-finest"? This should work. Or maybe "global-late" for evolution, since finest grid are evolved last? Yes, I think this should work.
As long as we have these modes, something like global-finest would be nice. I don't think it's a good idea to rely on the finest level being evolved last. That may be true by default, but not necessarily.
Frank
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
Hello all,
As long as we have these modes, something like global-finest would be nice. I don't think it's a good idea to rely on the finest level being evolved last. That may be true by default, but not necessarily.
We have GLOBAL, GLOBAL-EARLY and GLOBAL-LAST. I am not sure that adding GLOBAL-FINEST is needed. During EVOL the finest level has to be last because that is the way the mesh refinement works. As far as I remember GLOBAL during EVOL is GLOBAL-LAST for exactly the reason that interpolation then can be done in GLOBAL mode since all refinement levels have advanced so one can interpolate in time. I would suggest using GLOBAL-LAST over GLOBAL though since it makes it explicit what one wants.
Yours, Roland
- -- My email is as private as my paper mail. I therefore support encrypting and signing email messages. Get my PGP key from http://keys.gnupg.net.
On Mon, Oct 07, 2013 at 11:31:25AM -0700, Roland Haas wrote:
During EVOL the finest level has to be last because that is the way the mesh refinement works.
What I am saying is that while this is true at the moment in Carpet this doesn't need to be the case. Different levels can be evolved in any order as long as prolongation and restriction are taken care of in the correct order (e.g. by delaying them).
With grid arrays, being outside of the usual grid structure, it's not even really clear when and how often to evolve them. "Finest grid" - is that the maximally allowed finest grid, or the currently finest grid? If it is the latter - how do we deal with adding or removing a level? If it is the former - would global-late work, or wouldn't it skip steps for levels that the usual GFs aren't enabled?
Frank
On 7 Oct 2013, at 19:34, Erik Schnetter schnetter@cct.lsu.edu wrote:
On 2013-10-07, at 12:47 , Ian Hinder ian.hinder@aei.mpg.de wrote:
On 28 Sep 2013, at 16:17, Erik Schnetter schnetter@gmail.com wrote:
Yes, this (i.e. the work-around of scheduling the RHS in level mode, and then only actually executing the routine if this corresponds to global mode) is how I made things work. I haven't tried this recently, though. Maybe a test case would be in order...
Rather straightforward changes to MoL should allow scheduling the RHS in global mode as well, which would then work elegantly. This would probably require to duplicate MoL's integration logic, so that grid functions are integrated in level mode while grid arrays are integrated in global mode.
Is the problem that the grid arrays are not available in global mode, but they are available in level mode? Or is there some other reason that global mode cannot be used? When I try to call the interpolator in MoL_CalcRHS in level mode, I get
Assertion failed: (not (not want_global_mode and num_time_derivs.AT(m) == 0 and need_time_interp.AT(m))), function interpolate_components, file /Users/ian/Cactus/arrangements/Carpet/CarpetInterp/src/interp.cc, line 1353.
This is with InterpNumTimelevels = 1 in the interpolator parameters. I dug around in interp.cc, and came across the interpolation parameter want_global_mode. Setting this to 1 seems to allow the interpolation to succeed, but I don't know if this will break in certain situations. What does this interpolator parameter do?
My problem with registration turned out to be that I was trying to register MoL for grid SCALARs, not grid ARRAYs. A vector of grid arrays seems to work.
Grid arrays are available in global mode (and also in level, singlemap, and local mode).
Global mode cannot be used because MoL will still try to evolve them in level mode, i.e. more often. One needs to ensure that grid arrays have a non-zero RHS exactly on the finest grid. Global mode may be equivalent to the coarsest grid, which would then evolve with the wrong step size.
What does "equivalent to the coarsest grid" mean? Global mode routines are not called only every coarse grid timestep; I think they are called every fine grid iteration, when that level exists. As Frank pointed out, this could be problematic.
How/when are grid array timelevels rotated? I suspect this happens every fine grid iteration. What happens if a new finest grid is activated during a simulation? Does the effective dt of the grid array integration now become smaller? This means that the timelevels will have a nonuniform dt. Can MoL handle this?
Maybe there should be a concept of a grid array "belonging" to a specific refinement level, in the sense that it will only be "updated" when that refinement level exists. This means its timelevels will be rotated only then, and functions that want to access it will only be allowed to access it at those times. This would also allow integrations on a time step larger than the finest grid; at the moment, if I want to integrate something on the coarse grid, I have to do so with a timestep corresponding to the finest grid, which is obviously not very efficient.
The concept of having grid arrays separate from the normal grid functions is good in a "spatial" sense, but since Cactus treats timelevels globally, we have to deal with the issue of time refinement for grid arrays. An alternative would be to have a separate time hierarchy notion for grid arrays. Each grid array would have its own "dt". This might not even be commensurate with the gf time refinement factors. You would need a new way of scheduling a function "at" the time of a specific array (or group of arrays) in global mode, and during this call, you would likely need time interpolation if you wanted to access grid functions. I think this would be the best way of handling time stepping for grid arrays, but I think it would need some fairly serious changes to Carpet, and maybe also to parts of Cactus. Maybe attaching grid arrays to specific refinement levels is easier, and likely sufficient.
Yes, one cannot interpolate in level mode if this requires data from coarser levels. This is just what level mode means -- access data from the current refinement level only.
Is there a mode "global-finest"? This should work. Or maybe "global-late" for evolution, since finest grid are evolved last? Yes, I think this should work.
Do these exist at the moment, or would they need to be added?
Grid scalars are just a funny name for 0-dimensional grid arrays. If MoL doesn't support this, then it would be trivial to add -- the only difference is that they are called "scalar" instead of "array".
I'm not sure why MoL needs to even know about the group type (gf, scalar or array). I thought that Cactus flesh API was designed so that you can ask for a pointer to the variable, find the number of dimensions and size of the array, and loop over it, all without knowing whether it is a grid function or an array? Why does MoL know the difference? Knowing the difference means it has a lot of duplicated code, which is bad.
What does the interpolator parameter want_global_mode do? Should I be using it?
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
Hello Ian,
I'm not sure why MoL needs to even know about the group type (gf, scalar or array). I thought that Cactus flesh API was designed so that you can ask for a pointer to the variable, find the number of dimensions and size of the array, and loop over it, all without knowing whether it is a grid function or an array? Why does MoL know the difference? Knowing the difference means it has a lot of duplicated code, which is bad.
The issue comes from the scratch levels. For grid functions (which are all the same size but can only be allocated by the driver, though if we only need them for scratch one might not have to do so come to think of it) MoL uses the Scratch grid function. For grid arrays MoL uses malloc to get a chunk of memory itself (this is in InitialCopy.c). This chunk can differ from grid array to grid array since they all have different sizes.
Yes, MoL contains a lot of duplicated code. I hope to remove a lot of it after the IMEX code has been added. Though if someone else wants to clean it up instead you are very welcome.
Yours, Roland
- -- My email is as private as my paper mail. I therefore support encrypting and signing email messages. Get my PGP key from http://keys.gnupg.net.
On 7 Oct 2013, at 19:34, Erik Schnetter schnetter@cct.lsu.edu wrote:
On 2013-10-07, at 12:47 , Ian Hinder ian.hinder@aei.mpg.de wrote:
On 28 Sep 2013, at 16:17, Erik Schnetter schnetter@gmail.com wrote:
Yes, this (i.e. the work-around of scheduling the RHS in level mode, and then only actually executing the routine if this corresponds to global mode) is how I made things work. I haven't tried this recently, though. Maybe a test case would be in order...
Rather straightforward changes to MoL should allow scheduling the RHS in global mode as well, which would then work elegantly. This would probably require to duplicate MoL's integration logic, so that grid functions are integrated in level mode while grid arrays are integrated in global mode.
Is the problem that the grid arrays are not available in global mode, but they are available in level mode? Or is there some other reason that global mode cannot be used? When I try to call the interpolator in MoL_CalcRHS in level mode, I get
Assertion failed: (not (not want_global_mode and num_time_derivs.AT(m) == 0 and need_time_interp.AT(m))), function interpolate_components, file /Users/ian/Cactus/arrangements/Carpet/CarpetInterp/src/interp.cc, line 1353.
This is with InterpNumTimelevels = 1 in the interpolator parameters. I dug around in interp.cc, and came across the interpolation parameter want_global_mode. Setting this to 1 seems to allow the interpolation to succeed, but I don't know if this will break in certain situations. What does this interpolator parameter do?
My problem with registration turned out to be that I was trying to register MoL for grid SCALARs, not grid ARRAYs. A vector of grid arrays seems to work.
Grid arrays are available in global mode (and also in level, singlemap, and local mode).
Global mode cannot be used because MoL will still try to evolve them in level mode, i.e. more often. One needs to ensure that grid arrays have a non-zero RHS exactly on the finest grid. Global mode may be equivalent to the coarsest grid, which would then evolve with the wrong step size.
Yes, one cannot interpolate in level mode if this requires data from coarser levels. This is just what level mode means -- access data from the current refinement level only.
I don't follow this logic. In global mode, one can call the interpolator, even though there is no access to any data. CarpetInterp doesn't contain any macros that I can see for changing mode, and yet it calls functions for accessing data and sending it between processes. So it looks like CarpetInterp doesn't care much about the concept of modes. If I can call it in global mode, naively I don't see why I can't also call it in level mode and access data from the other levels. What is there that prevents this?
Is there a mode "global-finest"? This should work. Or maybe "global-late" for evolution, since finest grid are evolved last? Yes, I think this should work.
In unigrid, level mode appears to work. I haven't tried with refinement yet.
I have also tried global-late, but this does not work. With RK4, MoL_CalcRHS should be called four times, with cctk_time = (t^n, t^{n+1/2}, t^{n+1/2}, t^{n+1}) respectively. However, the first time it is called, it has cctk_time = 0, which is wrong. Calls to the interpolator lead to NaNs, which I also don't understand.
If I call it in level mode instead, then cctk_time is correct and the interpolator returns sensible-looking values. MoL_SetTime is called in level mode; is this related to the problem?
On 2013-10-16, at 9:09 , Ian Hinder ian.hinder@aei.mpg.de wrote:
On 7 Oct 2013, at 19:34, Erik Schnetter schnetter@cct.lsu.edu wrote:
On 2013-10-07, at 12:47 , Ian Hinder ian.hinder@aei.mpg.de wrote:
On 28 Sep 2013, at 16:17, Erik Schnetter schnetter@gmail.com wrote:
Yes, this (i.e. the work-around of scheduling the RHS in level mode, and then only actually executing the routine if this corresponds to global mode) is how I made things work. I haven't tried this recently, though. Maybe a test case would be in order...
Rather straightforward changes to MoL should allow scheduling the RHS in global mode as well, which would then work elegantly. This would probably require to duplicate MoL's integration logic, so that grid functions are integrated in level mode while grid arrays are integrated in global mode.
Is the problem that the grid arrays are not available in global mode, but they are available in level mode? Or is there some other reason that global mode cannot be used? When I try to call the interpolator in MoL_CalcRHS in level mode, I get
Assertion failed: (not (not want_global_mode and num_time_derivs.AT(m) == 0 and need_time_interp.AT(m))), function interpolate_components, file /Users/ian/Cactus/arrangements/Carpet/CarpetInterp/src/interp.cc, line 1353.
This is with InterpNumTimelevels = 1 in the interpolator parameters. I dug around in interp.cc, and came across the interpolation parameter want_global_mode. Setting this to 1 seems to allow the interpolation to succeed, but I don't know if this will break in certain situations. What does this interpolator parameter do?
My problem with registration turned out to be that I was trying to register MoL for grid SCALARs, not grid ARRAYs. A vector of grid arrays seems to work.
Grid arrays are available in global mode (and also in level, singlemap, and local mode).
Global mode cannot be used because MoL will still try to evolve them in level mode, i.e. more often. One needs to ensure that grid arrays have a non-zero RHS exactly on the finest grid. Global mode may be equivalent to the coarsest grid, which would then evolve with the wrong step size.
Yes, one cannot interpolate in level mode if this requires data from coarser levels. This is just what level mode means -- access data from the current refinement level only.
I don't follow this logic. In global mode, one can call the interpolator, even though there is no access to any data. CarpetInterp doesn't contain any macros that I can see for changing mode, and yet it calls functions for accessing data and sending it between processes. So it looks like CarpetInterp doesn't care much about the concept of modes. If I can call it in global mode, naively I don't see why I can't also call it in level mode and access data from the other levels. What is there that prevents this?
CarpetInterp certainly could interpolate from all levels; the question is what the user most likely wants, and to catch errors.
In level mode, the user most likely wants to operate on (and access) the current level only. Instead of interpolating from multiple levels, CarpetInterp assumes that it should interpolate only from the current refinement level. We later decided that interpolating from a single level wasn't a useful idea, and thus disabled it to catch more errors.
Global mode doesn't mean "cannot access grid functions"; it means "access all refinement levels of grid functions", as in "don't restrict yourself to a single level". As a by-product, one cannot define "pointer to THE grid function" since "THE" doesn't make sense. Thus pointers to individual grid functions are set to null in scheduled routines.
Yes, you probably want to use this interpolator parameter.
Is there a mode "global-finest"? This should work. Or maybe "global-late" for evolution, since finest grid are evolved last? Yes, I think this should work.
In unigrid, level mode appears to work. I haven't tried with refinement yet.
I have also tried global-late, but this does not work. With RK4, MoL_CalcRHS should be called four times, with cctk_time = (t^n, t^{n+1/2}, t^{n+1/2}, t^{n+1}) respectively. However, the first time it is called, it has cctk_time = 0, which is wrong. Calls to the interpolator lead to NaNs, which I also don't understand.
If I call it in level mode instead, then cctk_time is correct and the interpolator returns sensible-looking values. MoL_SetTime is called in level mode; is this related to the problem?
Likely yes. MoL modifies the current level's cctk_time, and not the global time. (Each level has its own time, and there is also a "global time" that corresponds to cctk_iteration * finest_delta_time.)
If MoL knew about global mode, and used it for grid arrays, then it would also need to modify the global time.
-erik
On 24 Sep 2013, at 11:34, Ian Hinder ian.hinder@aei.mpg.de wrote:
On 23 Sep 2013, at 18:06, Peter Diener diener@cct.lsu.edu wrote:
Hi,
No, you're right. The examples I mentioned are using either uni-grid or multi-patch without AMR. I didn't think of the complication of adding AMR to the mix.
So unless anything has changed, the latest information I have is from Erik's email from September 2007:
On 28 Sep 2007, at 18:10, Erik Schnetter schnetter@cct.lsu.edu wrote:
I am now evolving geodesics with the following work-around:
I use MoL. Each time the RHS is calculated, I check whether the current refinement level is the finest level. If so, the RHS is calculated, otherwise the RHS is set to zero.
The current refinement level is accessed by
interface.ccl: USES INCLUDE: carpet.hh
#include <carpet.hh>
if (Carpet::reflevel == Carpet::reflevels - 1) // do calculation
This works only from C++. The RHS calculating routine has to be scheduled in level mode.
I am using this method, and I have convinced myself that it is correct. The following facts are pertinent:
• Carpet loops through all iterations (of the finest possible grid), and within each iteration loops through all refinement levels from coarse to fine • For each iteration, Carpet only cycles the timelevels of grid arrays for the first (coarsest) refinement level that exists at that iteration. This means the timelevels are only cycled once, which is good. • MoL copies grid arrays from the past timelevel to the current timelevel on every refinement level
These facts combine to mean that no matter what RHS is used on the coarse grids, the correct past timelevel data is used to set the initial data for the iteration for MoL on the fine grid. I believe that you can compute the RHSs on all levels, without the check for the current level being the finest, and you still get the answer that was computed for the fine grid. This is because the values computed on the coarser grids are always overwritten (by MoL_InitialCopy) when the finer grid is evolved on that iteration. You could even set the RHSs to poison on the coarser grids (assuming that MoL can evolve poison without becoming slow or crashing). The latter is probably equivalent to scheduling the RHS routine in global-late, as I think Erik suggested in another post, as global-late routines are evolved only on the fine grid.
I suspect that this scheme will work if levels are added or removed during the simulation, but I haven't thought much about that case.
All the above ignores any potential interactions with the interpolator. I am currently calling the RHS routine in level mode, and setting the interpolator option "want_global_mode". There is a comment in CarpetInterp which says that this causes the interpolation to happen in global mode, i.e. with data from all timelevels. This seems to work; interpolation points not on the current level seem to give reasonable answers.
If I interpolate a variable within MoL_CalcRHS, I want the value at the current (MoL-set) value of cctk_time. This would be t_n, t_{n+1/2} or t_{n+1} for RK4. It looks like CarpetInterp uses the value of cctk_time to determine the time at which to interpolate. Interpolation performed on a coarse grid cannot give correct values from the fine grid if the interpolation point is in the fine grid, since the fine grid won't have been evolved to the required time yet. What happens in that situation? If Carpet aborts with an "extrapolation in time" error, then this means we have no choice but to evolve on the finest grid, which could be very inefficient.
I think we want to perform the interpolation in global-late (i.e. after all levels have been evolved), but not necessarily every iteration. Maybe we can fake this by setting the RHS to 0 on the iterations where we don't want to integrate, and setting it to rhs * dT/dt on the iterations where we do want to integrate, so that MoL will effectively use the correspondingly larger step size.
So the summary would be:
• Schedule the RHS routine in global-late, so it is only run on the finest level (and hence once all the others have been evolved). This probably causes MoL to integrate poison while on the other levels, but this should not make it into the output as the finest level will override the result. • Use the interpolator parameter "want_global_mode", as this allows the interpolation to use data from all refinement levels. • If you want to evolve less frequently than every fine-grid iteration, add the required logic to your RHS routine, and multiply the RHS by the corresponding factor.
Opinions?
On Nov 4, 2013, at 7:17 , Ian Hinder ian.hinder@aei.mpg.de wrote:
On 24 Sep 2013, at 11:34, Ian Hinder ian.hinder@aei.mpg.de wrote:
On 23 Sep 2013, at 18:06, Peter Diener diener@cct.lsu.edu wrote:
Hi,
No, you're right. The examples I mentioned are using either uni-grid or multi-patch without AMR. I didn't think of the complication of adding AMR to the mix.
So unless anything has changed, the latest information I have is from Erik's email from September 2007:
On 28 Sep 2007, at 18:10, Erik Schnetter schnetter@cct.lsu.edu wrote:
I am now evolving geodesics with the following work-around:
I use MoL. Each time the RHS is calculated, I check whether the current refinement level is the finest level. If so, the RHS is calculated, otherwise the RHS is set to zero.
The current refinement level is accessed by
interface.ccl: USES INCLUDE: carpet.hh
#include <carpet.hh>
if (Carpet::reflevel == Carpet::reflevels - 1) // do calculation
This works only from C++. The RHS calculating routine has to be scheduled in level mode.
I am using this method, and I have convinced myself that it is correct. The following facts are pertinent:
• Carpet loops through all iterations (of the finest possible grid), and within each iteration loops through all refinement levels from coarse to fine • For each iteration, Carpet only cycles the timelevels of grid arrays for the first (coarsest) refinement level that exists at that iteration. This means the timelevels are only cycled once, which is good. • MoL copies grid arrays from the past timelevel to the current timelevel on every refinement level
These facts combine to mean that no matter what RHS is used on the coarse grids, the correct past timelevel data is used to set the initial data for the iteration for MoL on the fine grid. I believe that you can compute the RHSs on all levels, without the check for the current level being the finest, and you still get the answer that was computed for the fine grid. This is because the values computed on the coarser grids are always overwritten (by MoL_InitialCopy) when the finer grid is evolved on that iteration. You could even set the RHSs to poison on the coarser grids (assuming that MoL can evolve poison without becoming slow or crashing). The latter is probably equivalent to scheduling the RHS routine in global-late, as I think Erik suggested in another post, as global-late routines are evolved only on the fine grid.
I suspect that this scheme will work if levels are added or removed during the simulation, but I haven't thought much about that case.
All the above ignores any potential interactions with the interpolator. I am currently calling the RHS routine in level mode, and setting the interpolator option "want_global_mode". There is a comment in CarpetInterp which says that this causes the interpolation to happen in global mode, i.e. with data from all timelevels. This seems to work; interpolation points not on the current level seem to give reasonable answers.
If I interpolate a variable within MoL_CalcRHS, I want the value at the current (MoL-set) value of cctk_time. This would be t_n, t_{n+1/2} or t_{n+1} for RK4. It looks like CarpetInterp uses the value of cctk_time to determine the time at which to interpolate. Interpolation performed on a coarse grid cannot give correct values from the fine grid if the interpolation point is in the fine grid, since the fine grid won't have been evolved to the required time yet. What happens in that situation? If Carpet aborts with an "extrapolation in time" error, then this means we have no choice but to evolve on the finest grid, which could be very inefficient.
I think we want to perform the interpolation in global-late (i.e. after all levels have been evolved), but not necessarily every iteration. Maybe we can fake this by setting the RHS to 0 on the iterations where we don't want to integrate, and setting it to rhs * dT/dt on the iterations where we do want to integrate, so that MoL will effectively use the correspondingly larger step size.
So the summary would be:
• Schedule the RHS routine in global-late, so it is only run on the finest level (and hence once all the others have been evolved). This probably causes MoL to integrate poison while on the other levels, but this should not make it into the output as the finest level will override the result. • Use the interpolator parameter "want_global_mode", as this allows the interpolation to use data from all refinement levels. • If you want to evolve less frequently than every fine-grid iteration, add the required logic to your RHS routine, and multiply the RHS by the corresponding factor.
Opinions?
Thank you for the analysis.
Note that MoL initializes the RHS to zero before calling CalcRHS. That is, the RHS will always be well-defined (no poison) by default.
-erik
users@lists.einsteintoolkit.org
-
Erik Schnetter -
Erik Schnetter -
Frank Loeffler -
Ian Hinder -
Peter Diener -
Roland Haas -
Yosef Zlochower