Incorrect totally trapped particle treatment when minimum in Bmag not at a grid point

[David Dickinson]

In invert_rhs_1 we consider a pitch angle to be the totally trapped particle if il == nlambda and there are trapped particles. This is already incorrect (at least in some sense) for the situation in which there are multiple trapped particles (see issue #235 for example).

If we have bmag which has > 1 local minima of different amplitudes then there should be multiple totally trapped particles. At least for Miller equilibria gridgen appears to only put a theta grid point at the theta=0 local minima (may not be a minima) whilst other local minima tend to be bracketed. In these bracketed cases the test to detect totally trapped particles, the a given pitch angle is forbidden at the left and right neighbouring grid points, does not flag these points as totally trapped. If these are in fact the global minima then the pitch angle which bounces here is il = nlambda yet if there is a (higher) local minima at theta=0 then there is a lower pitch index which is actually detected as a totally trapped particle but which is not considered totally trapped in invert_rhs_1.

The above is not very clearly phrased, but to summarise: In equilibria with multiple local minima we have incorrect treatment in invert_rhs_1 such that:

1. Trapped pitch angles which are totally trapped at some location but can travel along the field line at other theta (i.e. which are totally trapped at a local but not global minima in bmag) are not flagged as totally trapped and are just treated as normal trapped particles everywhere. We therefore rely on the regular trapped algorithm to ensure the totally trapped locations are dealt with correctly (i.e. that g(sigma_+) = g(sigma_-) at the totally trapped locations).
2. The largest pitch angle may not actually be considered totally trapped (as reflected in ittp) when the global minimum in bmag is not resolved by the theta grid (i.e. when there are neighbouring theta grid points bracketing the global minimum such that bmag is constant between two neighbouring grid points) but it is considered totally trapped in invert_rhs_1.

PR #767 attempts to improve the handling of totally trapped particles for this and issue #235.

[David Dickinson]

Here’s an example geometry which has multiple local minima with the global minimum on the grid at theta = 0

&theta_grid_parameters                                                                                                                                                                        
 ntheta = 128                                                                                                                                                                                 
 nperiod = 1                                                                                                                                                                                  
 shift = 0.15                                                                                                                                                                                 
 rhoc = 0.95                                                                                                                                                                                  
 akappa = 3.0                                                                                                                                                                                 
 akappri = -2.5                                                                                                                                                                               
 tri = 0.5                                                                                                                                                                                    
 tripri = 0.0                                                                                                                                                                                 
/                                                                                                                                                                                             
&theta_grid_knobs                                                                                                                                                                             
equilibrium_option = "eik"                                                                                                                                                                    
 /                                                                                                                                                                                            
 &theta_grid_eik_knobs                                                                                                                                                                        
 local_eq = T                                                                                                                                                                                 
bishop = 4                                                                                                                                                                                    
beta_prime_input = -1.0                                                                                                                                                                       
writelots = T                                                                                                                                                                                 
s_hat_input = 0.8                                                                                                                                                                             
ntheta_geometry = 4096                                                                                                                                                                        
equal_arc = F                                                                                                                                                                                 
 /        

‌

[David Dickinson]

Here’s an example geometry which has multiple local minima in bmag but with the global minimum not resolved by the theta grid.

&theta_grid_parameters                                                                                                                                                                        
 ntheta = 128                                                                                                                                                                                 
 nperiod = 1                                                                                                                                                                                  
 shift = -0.3                                                                                                                                                                                 
 rhoc = 0.95                                                                                                                                                                                  
 akappa = 3.0                                                                                                                                                                                 
 akappri = 2.5                                                                                                                                                                                
 tri = -0.5                                                                                                                                                                                   
 tripri = 1.0                                                                                                                                                                                 
/                                                                                                                                                                                             
&theta_grid_knobs                                                                                                                                                                             
equilibrium_option = "eik"                                                                                                                                                                    
 /                                                                                                                                                                                            
 &theta_grid_eik_knobs                                                                                                                                                                        
 local_eq = T                                                                                                                                                                                 
bishop = 4                                                                                                                                                                                    
beta_prime_input = -1.0                                                                                                                                                                       
writelots = T                                                                                                                                                                                 
s_hat_input = 0.8                                                                                                                                                                             
ntheta_geometry = 4096                                                                                                                                                                        
equal_arc = F                                                                                                                                                                                 
 /                   

‌

[David Dickinson]

Here’s a comparison of the two cases (first in orange, second in blue):

‌

[David Dickinson]

Running a simple test case with next and PR #767 for the two equilibria above show very small differences suggesting the generalised treatment is not significant here (although it should be noted the solution is not particularly well behaved). One thing to note is that the improved treatment still relies on ittp so does not change the treatment around bracketed minima but only changes the behaviour for resolved totally trapped points with bmag > bmin. Here that means we’d only expect an effect in the blue equilibria.

In addition to PR #767 we should check that the regular trapped algorithm correctly deals with bracketed minima such as those for theta/=0 in the orange case above.