[Carpet] Prolongation order

[Erik Schnetter]

On Apr 2, 2008, at 17:23:06, Ian Hinder wrote:
> Hi,
>
> I have a question about prolongation operators in Carpet.  It's very
> possible that I have some misunderstandings in what follows.  For the
> operator given by prolongation_order_space=5, I assume that the  
> relevant
> source file is CarpetLib/src/prolongate_3d_o5_rf2.cc.

This is correct.

> From this file, which I do not fully understand, it looks to me like  
> the
> interpolation uses the points i-2, i-1, i, i+1, i+2, i+3.  In other
> words, when interpolating to any given point, there is a preferred
> direction.

Yes, these are the points which are used.  However, the indexing is  
different, since the different levels have different grid spacing, and  
the point corresponding to "i" does not exist on the coarser grid.   
The stencil is supposed to consist of 6 point, 3 strictly below and 3  
strictly above the target point.

>  Shouldn't this mean that the interpolation is no longer
> invariant when (for example) a reflection symmetry thorn is used?  A
> polynomial interpolation which uses 6 points has a fitting  
> polynomial of
> degree 5 (this gives the order in the parameter), and an interpolation
> error of order 6.  Is this correct so far?

This is what I call a 5th order operator, and this is what is supposed  
to be used.  I want to use stencils with even number of points to keep  
it symmetric.

> Now, if I want to perform an interpolation, I assume that the nearest
> grid point to the coordinate that I want to interpolate to will be
> labeled i, and I will need at least three points in a neighbourhood  
> of i
> to be available to perform the interpolation.  In CarpetLib/src/dh.cc,
> we have
>
> prolongation_stencil_size () = ... prolongation_order_space / 2
>
> which I understand will give 2, as C rounds integer arithmetic  
> down.  I
> would have expected to have 3, as the required number of points.  Is
> this not the intent of this parameter?

Yes, and no.  The calculation is somewhat complex; however, this  
function is indeed badly named.

A fine grid near a refinement boundary may end with a grid point which  
does not exist on the next coarser grid.  In order to calculate the  
source points for a prolongation stencil, the fine grid is first  
"expanded_for" the next coarser grid, i.e., it is "rounded outwards",  
so to say.  This may add one additional grid point.  Alternatively, if  
the fine grid ends with a grid point which does exist on the next  
coarser grid, then this grid point does not need interpolation (it is  
copied instead).

This means that one point of the interpolation stencil is always there  
"for free", and therefore prolongation_stencil_size is smaller by  
one.  This should arguably be handled by increasing  
prolongation_stencil_size and subtracting one explicitly.  However,  
the current expression is also correct for cell-centred grids where  
the orders are even numbers.

Note that prolongation_stencil_size defines only the size of the  
source region for a prolongation operation.  Each prolongation  
operator checks the sizes again independently to ensure that  
inconsistencies are caught.  The number of required ghost zones on the  
coarser grids is (prolongation_order_space+1)/2.  This expression is  
also correct for cell-centred grids.

> Finally, since the number of ghost zones required is 3, we could
> implement 6th order prolongation (with a 7th order error) for a small
> computational cost, without requiring any more ghost zones, and  
> without
> introducing a preferred direction.

The stencil which is called 5th order in Carpet requires 3 ghost zones.

-erik

-- 
Erik Schnetter <schnetter at cct.lsu.edu>   http://www.cct.lsu.edu/~eschnett/

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: 194 bytes
Desc: This is a digitally signed message part
Url : /archives/developers/attachments/20080402/02d4ac81/attachment.pgp