ideal_ball error dependant on no. of beta_prime/shat values

[Former user]

ideal_ball will crash with the following error if the no. of beta_prime/shat values is increased too high

 FITPACK: curv1 error: x-values not increasing
application called MPI_Abort(MPI_COMM_WORLD, 13) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=13
:
system msg for write_line failure : Bad file descriptor

‌

I’ve attached a simple case (GA-STD like) where this will fail with 10 shat and 10 beta_prime values. The error seems to occurs with a large shift (in this case -0.3) and will disappear if shift is dropped. The actual values on the grid don’t seem to matter, just the number of gridpoints. This case runs fine with nshat=10 and nbeta_prime=6 and vice versa.

‌

The actual error arises in the geometry.f90 file when GS2 tries to interpolate the theta_prime into an equally space grid (in interp_cspl subroutine). When I plot the theta_prime that goes into the function I get the following and there’s no discernible difference between them and they both are continuously increasing.

I’m not entirely sure if theta_prime should look like this anyway. Should it be a smooth increasing value unless I’ve misunderstood its use here.

‌

[David Dickinson]

My guess is that something isn’t being cleaned up/recalculated between repeat calls to eikcoefs. I know @Peter Hill was looking at similar issues recently.

[Former user]

So the theta grid isn’t being initialised correctly. Below I’ve plotted the theta grid for each call to eikcoeff. Each iteration seem to have the points clustered towards a multiple of pi

‌

[Former user]

Setting equal_arc to .false. seems to fix this issue.

With equal_arc = .true. each iteration calling eikcoeff seems to redefine the theta grid to try and be equally spaced in theta_prime. It should produce the same theta grid after the first iteration so there’s likely an issue with the implementation of equal_arc

[David Dickinson]

It looks like the theta grid is retaining the equal_arc derived values from previously calls to eikcoefs but the code around equal arc assumes that theta is a uniformly spaced grid of theta values (or at least the comments say this is the case).

‌

[David Dickinson]

OK I think I’ve found the issue. This input file uses the Miller geometry type, for which the theta grid seems to actually be calculated from a call to init_theta_grid unlike all other eik geometries for which the theta grid is calculated in eikcoefs as well. Because the s-alpha scan is essentially just calling eikcoefs many times for other eik types each call resets the theta grid but for Miller we’re left with the theta values that we had at the end of the previous call rather than a fresh theta grid. A hacky fix for now is to call init/finish_theta_grid at the start/end of check_stability in the ballstab module. Something like

‌

‌

  subroutine check_stability(shat_in, beta_prime_in, ib, iunstable, restore)
    !Note iunstable>0 means unstable - Could use a logical but size may be interesting                                                                                                        
    use geometry, only: dbetadrho, eikcoefs
    use theta_grid, only: init_theta_grid, finish_theta_grid
    implicit none
    real, intent(in) :: shat_in, beta_prime_in
    integer, intent(in) :: ib
    integer, intent(out) :: iunstable
    logical, intent(in), optional :: restore
    real :: shat_bak, beta_prime_bak
    call init_theta_grid
    !Store initial state                                                                                                                                                                      
    call get_shat(shat_bak)
    call get_beta_prime(beta_prime_bak)

    !Setup geometry                                                                                                                                                                           
    !NOTE: As how we can influence beta_prime (and shat) varies with bishop                                                                                                                   
    !we actually end up setting different parameters in different cases. The                                                                                                                  
    !physical parameter dbetadrho should really be the relevant parameter.                                                                                                                    
    call set_beta_prime(beta_prime_in)
    call set_shat(shat_in)
    call eikcoefs !Recalculate geometry terms                                                                                                                                                 

    !Here we store the value of dbetadrho, this should be bishop value independent                                                                                                            
    dbdrho_arr(ib)=dbetadrho

    !Test stability                                                                                                                                                                           
    iunstable=is_unstable()

    !Restore geometry if requested                                                                                                                                                            
    if(present(restore))then
       if(restore) then
          call set_beta_prime(beta_prime_bak)
          call set_shat(shat_bak)
          call eikcoefs
       endif
    else
       call finish_theta_grid
    endif

  end subroutine check_stability

‌

This probably isn’t a good long term fix (maybe have to move the init/finish theta call to the 2D loop rather than here).

‌

Another alternative is to make sure that the Miller theta grid is actually initialised in eikcoefs

[Former user]

If there are other scenarios in which eikcoeff is called multiple times then its probably best to put the init_theta_grid inside eikcoeff. Seems like a safer bet

‌

[Peter Hill]

Think I have a cleaner solution, will make a PR tomorrow hopefully

[David Dickinson]

Yes that’s the best fix but will require a bit more work to get it in.

[Peter Hill]

Fix initial theta grid size potentially unset in call to eikcoefs

Fixes #111

For Miller/local and EFIT equilibria only, it was necessary to call geometry::init_theta before the call to eikcoefs, while for the remaining numerical equilibria, a call in eikcoefs to tdef was necessary after reading the equilibrium file.

Additionally, tdef and init_theta were essentially identical, with the exception that tdef forced ntheta to be even.

This commit will:

• Replace tdef and init_theta with init_uniform_theta_grid

• Always call init_uniform_theta_grid inside eikcoefs in the appropriate places

• Make init_uniform_theta_grid private; it is now always called appropriately, so no other calls are needed

• Replace duplicated calls to tdef with single call to init_uniform_theta_grid

• Make the ntheta_returned out-parameter to ekicoefs no longer optional, and make it intent(inout); callers now always get the updated value

• Add docstrings to various touched variables/procedures, making sure they are all consistent

→ <<cset ad84d40722a3>>

[David Dickinson]

Fix initial theta grid size potentially unset in call to eikcoefs

Fixes #111

For Miller/local and EFIT equilibria only, it was necessary to call geometry::init_theta before the call to eikcoefs, while for the remaining numerical equilibria, a call in eikcoefs to tdef was necessary after reading the equilibrium file.

Additionally, tdef and init_theta were essentially identical, with the exception that tdef forced ntheta to be even.

This commit will:

• Replace tdef and init_theta with init_uniform_theta_grid

• Always call init_uniform_theta_grid inside eikcoefs in the appropriate places

• Make init_uniform_theta_grid private; it is now always called appropriately, so no other calls are needed

• Replace duplicated calls to tdef with single call to init_uniform_theta_grid

• Make the ntheta_returned out-parameter to ekicoefs no longer optional, and make it intent(inout); callers now always get the updated value

• Add docstrings to various touched variables/procedures, making sure they are all consistent

→ <<cset 866273646ab5>>