Differences between flux tube and ballooning space simulations

[David Dickinson]

One might expect linear flux tube and ballooning space simulations, appropriately set up, to be entirely equivalent.

As demonstrated in issue #248, this is not always the case in all areas. This issue is intended to hold discussion of differences etc.

‌

[David Dickinson]

Gridgen can produce slightly different grids between flux tube with nperiod = 1 and nperiod > 1.

[David Dickinson]

The wfbbc_option = 'trapped' approach only treats the end of a cell as a wfb bounce point. This means that a linked flux tube has bounce points per cell whilst a ballooning space run would only have bounce points at the far ends of the supercell.

‌

[David Dickinson]

Looking at a simple CBC like test case I have been comparing the distribution function. Whilst there’s mostly good agreement the WFB shows more obvious differences. Here’s the magnitude of the distribution function on a log scale for ions with the lowest energy travelling from left to right. Red dashed is the unlinked nperiod = 2 case whilst the other colours are from a linked flux tube case with nperiod =1

‌

One can see clear differences near the incoming boundary. Here’s the same comparison for the passing pitch angle just before the wfb

The agreement is much better (although it still has the sharp oscillation rather than the smoother solution.). These runs are with nonad_zero = T. Here’s the same comparison with nonad_zero = F for the wfb

and the passing pitch angle just before the wfb

These runs use wfbbc_option = 'passing'

[David Dickinson]

It seems the difference for the wfb between the nperiod = 2 and nperiod = 1 cases above are sensitive to the choice of the boundary option. For the nperiod = 2 case I run with nx = 1 and zero boundaries. This should be equivalent to linked boundaries as nx = 1 yet here is the nonad_zero = F wfb comparison as above but using linked for both cases.

It appears that the difference is due to the is_self_periodic flag in invert_rhs_1. This is set through

       if (boundary_option_switch == boundary_option_linked) then
          is_self_periodic = (speriod_flag .and. is_passing) & ! Passing particle with ky=0                                                                                                   
               .or. (speriod_flag .and. is_wfb .and. .not. trapped_wfb) & ! Non-trapped wfb treatment of wfb with ky = 0                                                                      
               .or. (is_wfb .and. .not. connections(iglo)%neighbor & ! Isolated (i.e. ky/=0 cell with no connections) wfb                                                                     
               .and. mixed_wfb )     ! treated as neither passing or trapped particle (i.e. unique)                                                                                           
       else
          is_self_periodic = (use_pass_homog .and. is_passing) & ! Passing particle with ky=0 or with self-periodic boundaries                                                                
               .or. (is_wfb .and. .not. trapped_wfb) ! Non-trapped treatment of wfb with any B.C.                                                                                             
       end if

For these runs this evaluates to true in the nperiod = 2, boundary_option = ‘zero’ run whilst it evaluates to false in the nperiod = 2, boundary_option = ‘linked’ case. For consistency I think the non-linked branch should be

         is_self_periodic = (use_pass_homog .and. is_passing) & ! Passing particle with ky=0 or with self-periodic boundaries                                                                
               .or. (is_wfb .and. mixed_wfb) ! Non-trapped treatment of wfb with any B.C.       

as it is only the mixed_wfb treatment which is supposed to be periodic on the supercell.

Of course we may note that the periodic treatment of the wfb appeared to avoid the sharp oscillations at the boundary. This is also avoided when using the “correct” trapped wfb treatment. This perhaps serves to demonstrate that it is the zero incoming boundary conditions which are (at least partially) responsible for the structure that we see near the boundary here.

It should be noted that the plots shown here are of gnew (i.e. output from save_for_restart = T) and not the non-adiabatic dfn h (i.e. the output of save_distfn = T).

[David Dickinson]

To help demonstrate that the oscillation may not be due to the relatively low theta extent used in the previous test case here is the distribution function for the passing pitch angle just before the wfb for a case with nx = 15 (or nperiod = 8) and nonad_zero = F

Here’s the same with nonad_zero = T

This persists when looking at the non-adiabatic part h (nonad_zero = T):

and indeed no similar structure is observed in the phi/apar structure in this linear run.

This may be a consequence of the initial conditions – here I’m using ginit_option = ‘default' which is known to not satisfy the initial conditions, however similar structure can still be observed in runs using ginit_option = 'random_sine':

‌

[David Dickinson]

One thought I’ve had previously about some of the known issues around the boundary condition enforcement is that enforcing a boundary condition exactly on the grid is known to cause problems in other numerical approaches. Here I’ve added in a pink ballooning space case which attempts to enforce gnew == 0 exactly half a theta grid point beyond the end of the grid. I do this by evaluating the parallel solve starting with gnew set the usual way at the end of the grid to find gnew one point in from the boundary, use these two values of gnew to form the central difference approximation of dg/dtheta at the grid centre and then extrapolate back half a grid point beyond the end of the grid. Setting this to zero one can then rearrange to find out the value of gnew on the boundary point which satisfies this. Explicitly the code looks like:

       ! Integrate along theta one grid point 
       ig = ntgrid - 1
       gtmp = -gnew(ig+1,2,iglo)*r(ig,2,iglo) + ainv(ig,2,iglo)*source(ig,2,iglo)
       ! We can approximate the derivative at ig + 1/2 as [gtmp - gnew(ntgrid,2,iglo)]/dtheta.
       ! We can the extrapolate to ig = ntgrid + 1/2 as g_{ig+1} = g_boundary + (dtheta/2) * dg/dtheta
       ! --> g_{ig+1/2} ~ g_boundary + [gtmp - gnew(ntgrid,2,iglo)] * (dtheta/2) / dtheta
       ! Now set g_{ig+1/2} = 0 and solve for g_boundary --> 
       ! g_boundary = 0.5 * (gnew(ntgrid,2,iglo) - gtmp)
       gnew(ntgrid, 2, iglo) = (gtmp - gnew(ntgrid, 2, iglo)) * 0.5

       ! inhomogeneous part: gnew                                                                                                                                                             
       ! r=ainv=0 if forbid(ig,il) or forbid(ig+1,il), so gnew=0 in forbidden                                                                                                                 
       ! region and at upper bounce point                                                                                                                                                     
       do ig = ntgrid-1, -ntgrid, -1
          gnew(ig,2,iglo) = -gnew(ig+1,2,iglo)*r(ig,2,iglo) + ainv(ig,2,iglo)*source(ig,2,iglo)
       end do

The new pink case appears to show much fewer (but not zero) grid scale oscillations. The above is with nonad_zero = start_from_previous_solution = F whilst the below sets both these flags to true and can be seen to suffer from a bit more grid scale in from the boundary.

It seems that nonad_zero = T contributes to this a bit but start_from_previous_solution is the main cause. Below is nonad_zero = F ; start_from_previous_solution = T

‌

Edit:

Replacing gnew(ntgrid, 2, iglo) = (gtmp - gnew(ntgrid, 2, iglo)) * 0.5 with gnew(ntgrid, 2, iglo) = (gtmp - gnew(ntgrid, 2, iglo)) * 0.5 * (1 + bkdiff(is)) seems to further suppress the boundary grid scale for the lowest energy ions but things worse for the higher energy ions. Perhaps this is indicating that we should weight the extrapolation differently.

[David Dickinson]

The previous plots were all for the lowest energy ion. Taking the most energetic ion (with off grid boundaries) we can still identify less pleasing structure:

Interestingly this is less visible in the electrons, although their structure also has some concerning features:

Here’s the structure for the most energetic ion for a trapped pitch angle, which also shows troubling structure.

‌