Prolongation/Restriction of scalars in foreach_neighbor/stencils

[Graham Garcia]

Hi all,

I am working on a function that can output VOF interface points with their normal values more "smooth" to each other by utilizing foreach_neighbor() to generate a stencil of interface values at each interface and smoothing accordingly. however, the working approach will calculate each interface point at each point in the foreach_neighbor, meaning that in simulations with lots of interfacial data, will be re-calculating the same thing over and over again. To prevent this I create a scalar that tracks the 'ID-location' of the point, and an attributed structure to the field which builds and computes the vof interface, and the smoothened vof interface for each scalar with an "ID", note all of these are marked at interficial locations. 

This method works well until using a foreach_neighbor, as the ID will get smoothened by Basilisk's default prolongation and restriction techniques if it exists on a separate level. So I apply the following methods to prevent the 'smoothening' of an ID

-For prolongation, i find the location of its point and add its ID there

-For restriction, I take the first real child value and set it to the parent.

 #if TREE
 //when prolongating we use our structure to find the location 

static void smoothvof_prolongation(Point point, scalar s){
     double val = s[];
     double *ppoint = s.smooth->points[(int)s[]];
     foreach_child(){
         //We need to refrence the structure to determine which cell will hold our lower ID
         if(insidebounds(ppoint,x,y,Delta)){
             s[] = val;
         }
         else{
             s[] = nodata;
         }
     }
 }

 static inline void smoothvof_restriction(Point point, scalar s){
     int pass = 0;
     double val = nodata;
     foreach_child(){//we pick the first child
         if(!pass && s[] != nodata){
             val = s[];
             pass++;
         }
     }
     s[] = val;
 }
 #endif

I then set my scalar & its prolongation  and restriction techniques

  scalar vofref[];

  vofref.refine = vofref.prolongation = smoothvof_prolongation;
  vofref.coarsen= vofref.restriction  = smoothvof_restriction;

we can then loop and set all values. 

and accessing them later using our ref value

foreach(){

  if(c > 1e-6 && c < 1-1e-6 && vofref[] != nodata){

      //do some allocation

      int ref = (int)vofref[];

    foreach_neighbor(){

      if(c > 1e-6 && c < 1-1e-6 && vofref[] != nodata){

        int newref = (int)vofref[];

      }

   }

  }

}

my value at 'ref' will always be correct, and have the ID it is meant to, but inside the foreach_neighbor and looking at deeper levels, I am seeing fractional values for my 'newref' which is wrong. 

This leads me to my questions:

-How does basilisk ensure a correct stencil inside a foreach() when dealing with different computation levels? does it use prolongation&restriction or something separate? Placing an output inside of my prolongation and restriction methods, I can see that they are never being utilized within/before a foreach loop, yet the stencil is still guaranteed and being generated somehow. 

-Furthermore does the restriction method only work when t=0.0? I place a restriction({vofref}) after setting my vofref IDs, But placing an output in the restriction method only shows up at t=0.000, and the other timesteps where the vofref is computed the same as before it never gets run at all. The scalar is created and cleaned within an event at each output timestep, so I feel like I should see restriction running at every timestep.

My current method does computation using MPI, however, all my issues arise from the incorrect prolongation/refinement on my scalar holding my IDs, so I doubt it is related to the errors.

Any help would be much appreciated!

[j.a.v...@gmail.com]

Hallo,

> How does basilisk ensure a correct stencil inside a foreach() when dealing with different computation levels?

It makes a V-Cycle, first using the a restriction (indeed not by calling the function "restriction()") from high to low levels, and then applying the prolongation from low to heigher levels.  For details see: the function "tree_boundary_level()" in /src/grid/tree-common.h

> Furthermore does the restriction method only work when t=0.0?

Maybe, see the function "tree_restriction()" in the aforementioned source file. If you want to enforce restriction for debugging, you should call "multigrid_restriction()".

When I dealt with a similar issue, I stored the coordinates/normals in a vector field.

Antoon

Op dinsdag 31 oktober 2023 om 21:05:30 UTC+1 schreef graha...@gmail.com:

[Graham Garcia]

Hi,

I see, if I'm correct the foreach itself doesn't enforce the tree_boundary_level() itself, as if I run boundary_level({vofref},max_level) then it can do the prolongation and restriction that I defined. vofref is defined and deleted inside a function, so I'm guessing that's why the prolongation and restriction were never showing up. 

Thanks!

Graham Garcia