ideal_ball error dependant on no. of beta_prime/shat values
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.
Comments (10)
-
-
Account Deleted reporter 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
-
Account Deleted reporter 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
-
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).
-
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 othereik
geometries for which the theta grid is calculated ineikcoefs
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 ofcheck_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
-
Account Deleted reporter 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
-
Think I have a cleaner solution, will make a PR tomorrow hopefully
-
Yes that’s the best fix but will require a bit more work to get it in.
-
- changed status to resolved
Fix initial theta grid size potentially unset in call to eikcoefs
Fixes
#111For Miller/local and EFIT equilibria only, it was necessary to call
geometry::init_theta
before the call toeikcoefs
, while for the remaining numerical equilibria, a call ineikcoefs
totdef
was necessary after reading the equilibrium file.Additionally,
tdef
andinit_theta
were essentially identical, with the exception thattdef
forcedntheta
to be even.This commit will:
-
Replace
tdef
andinit_theta
withinit_uniform_theta_grid
-
Always call
init_uniform_theta_grid
insideeikcoefs
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 toinit_uniform_theta_grid
-
Make the
ntheta_returned
out-parameter toekicoefs
no longer optional, and make itintent(inout)
; callers now always get the updated value -
Add docstrings to various touched variables/procedures, making sure they are all consistent
→ <<cset ad84d40722a3>>
-
Fix initial theta grid size potentially unset in call to eikcoefs
Fixes
#111For Miller/local and EFIT equilibria only, it was necessary to call
geometry::init_theta
before the call toeikcoefs
, while for the remaining numerical equilibria, a call ineikcoefs
totdef
was necessary after reading the equilibrium file.Additionally,
tdef
andinit_theta
were essentially identical, with the exception thattdef
forcedntheta
to be even.This commit will:
-
Replace
tdef
andinit_theta
withinit_uniform_theta_grid
-
Always call
init_uniform_theta_grid
insideeikcoefs
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 toinit_uniform_theta_grid
-
Make the
ntheta_returned
out-parameter toekicoefs
no longer optional, and make itintent(inout)
; callers now always get the updated value -
Add docstrings to various touched variables/procedures, making sure they are all consistent
→ <<cset 866273646ab5>>
-
- Log in to comment
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.