[Carpet] question on grid structures and interpolation

[Erik Schnetter]

Hi Yosef,

if you use MoL, you should probably use either local or level mode to  
calculate the RHS.  This is because MoL will evolve the variable in  
level mode.

You can schedule a routine in level mode which sets the RHS to zero,  
then schedule another routine in global mode after the first routine  
which calculates the RHS correctly.  In this way you can fake a  
continuous time evolution in spite of using a time integrator which  
is discontinuous in time (i.e., uses the Berger-Oliger scheme).

There is no good way to find out whether a level contains a point,  
except to look at Carpet's data structures.  Be aware that this can  
change during regridding, and regridding can happen between the  
coarse and fine grid time steps.

The complication here is conceptual.  One has to be very careful to  
define how the fine grid should be evolved in the presence of the  
discontinuous Berger-Oliger time evolution, with a potentially  
changing grid structure underneath.

I often define several time levels manually (i.e., not using the  
TIMELEVELS specifier in the interface), and I also store explicitly  
the time for each time level.  This allows changes in the time step  
size, which can happen if the grid structure changes.  I then  
implement my own time stepping algorithm, usually simple and second- 
order accurate, using 3 time levels altogether.  This is all called  
in global mode.  Not tying it to MoL or level mode ensures that  
everything remains consistent.

Something similar could be implemented in MoL, if MoL were to evolve  
grid arrays in global mode.  Note that the coupling between evolved  
grid functions and grid arrays is then still only first order  
accurate.  I believe that anything better would be highly non-trivial.

-erik

On Aug 30, 2007, at 09:31:08, Yosef Zlochower wrote:

> Hi,
>
>  I am trying to evolve curves using MoL. I have cactus arrays
> loc_x, loc_y, loc_z, and  dot_loc_x, dot_loc_y, dot_loc_z. I then
> have dot_loc_[xyz] = INTERPOLATION(ADMBASE_VARS).
> MoL will then handle the time evolution. Do I need to schedule this
> routine in local mode for MoL to work correctly? Is there a way of  
> finding
> if the current level actually contains some point p?
>
>
>
> Yosef
>
> Erik Schnetter wrote:
>> On Aug 30, 2007, at 08:40:19, Yosef Zlochower wrote:
>>
>>> Hi,
>>>
>>>   I have a routine, which is scheduled in local mode, that
>>> interpolated GFs
>>> at specific points. I wanted to set up the routine so that it  
>>> would only
>>> do the interpolation if the current component of the current level
>>> contains
>>> the interpolation point. I tried to get the bounds of the current
>>> component
>>> using the code
>>>
>>> lx  = CCTK_ORIGIN_SPACE(0);
>>> ux = lx + (cctk_gsh[0] -1)*CCTK_DELTA_SPACE(0);
>>> etc.
>>>
>>> My grid is set up so that I have n refinement levels about two  
>>> separate
>>> locations
>>> (using the CarpetRegrid2 thorn). I found, however, that the range
>>> (lx,ux) seems to
>>> contain both centers for all refinement levels. Does this mean  
>>> that I am
>>> actually
>>> using much bigger grids than I thought, or does the above  
>>> calculation
>>> only give
>>> me results for a `virtual' grid which contains all components of the
>>> real grids at
>>> a particular refinement level?
>>
>> Hi Yosef,
>>
>> the information stored in the cGH can only represent a cuboid
>> (rectangular) box.  It therefore represents a box which is guaranteed
>> to enclose all components of the current level, but not all grid
>> points will actually exist as components.  Yes, this is a  
>> "virtual" box.
>>
>>> In addition, would be the result if I attempt to interpolate a data
>>> point from another
>>> component of the same refinement level?
>>
>> Performing interpolations in local mode can lead to strange results.
>> If there are several components per processor, then the  
>> interpolations
>> will occur multiple times.  They are mostly allowed for backward
>> compatibility with code which doesn't know about modes -- you should
>> try to use either global or level mode instead.  Level mode
>> interpolates from the current level, which seems to fit your problem.
>>
>> Interpolations are always global, so you can always get data from
>> other processors (other components). All processors need to call the
>> interpolation routines at the same time.
>>
>> If you really only want to interpolate if the current (this
>> processor's) component of the current level contains the  
>> interpolation
>> point, then I suggest to use a local interpolation instead.  A local
>> interpolation does not know about symmetry boundaries, but it is
>> faster since it does not involve any communication, and it does not
>> have to be called in unison on all processors.  This is done by
>> calling CCTK_InterpLocalUniform, which expects instead of grid
>> function indices a description of the coordinate system and pointers
>> to the arrays which you want to interpolate.
>>
>> In my language, a component is a cuboid 3D piece of data residing  
>> on a
>> single processor.  If you have two "refined regions", then these will
>> still be broken up into multiple components.  Note that these won't
>> necessarily be two separate sets of components; if you have e.g. 2
>> regions and 2 processors, then each region will become one component.
>>
>> -erik
>>
>> --Erik Schnetter <schnetter at cct.lsu.edu>
>>
>> My email is as private as my paper mail.  I therefore support  
>> encrypting
>> and signing email messages.  Get my PGP key from www.keyserver.net.
>>
>>
>>
>> --------------------------------------------------------------------- 
>> ---
>>
>> _______________________________________________
>> developers mailing list
>> developers at lists.carpetcode.org
>> http://lists.carpetcode.org/listinfo/developers
>>
>
> _______________________________________________
> developers mailing list
> developers at lists.carpetcode.org
> http://lists.carpetcode.org/listinfo/developers
>


-- 
Erik Schnetter <schnetter at cct.lsu.edu>

My email is as private as my paper mail.  I therefore support encrypting
and signing email messages.  Get my PGP key from www.keyserver.net.



-------------- next part --------------
A non-text attachment was scrubbed...
Name: PGP.sig
Type: application/pgp-signature
Size: 186 bytes
Desc: This is a digitally signed message part
Url : /archives/developers/attachments/20070901/18aaa713/attachment.pgp