Skip to content

guess_axis scans only half the toroidal planes when lasym = T #21

Description

@CharlesCNorton

guess_axis computes ivmax for the two symmetry cases,

if (.not.lasym) then
  ivmax = nzeta/2+1
else
  ivmax = nzeta
end if

and then never uses it as the loop bound: the plane loop is planes: DO iv = 1, nzeta/2+1 (guess_axis.f90:72). For lasym = F the remaining planes are filled after the loop by the stellarator-symmetry mirror

do iv = nzeta/2+2, nzeta
    rcom(iv) = rcom(nzeta+1-(iv-1))
    zcom(iv) =-zcom(nzeta+1-(iv-1))
end do

which is correctly skipped for lasym = T. Nothing else writes rcom and zcom, so for lasym = T the entries nzeta/2+2 .. nzeta reach the axis Fourier transform

raxis_cc(n) = dzeta*SUM(cosnv(:,n)*rcom(:))/nscale(n)

never having been assigned. Both arrays are allocatable, so no bounds or initialization check catches it.

PARVMEC guess_axis.f:132 runs its loop over all nzeta planes and uses the mirror as an in-loop shortcut for .not.lasym instead:

planes: DO iv = 1, nzeta
   IF (.not.lasym .and. iv .gt. nzeta/2 + 1) THEN
      rcom(iv) = rcom(nzeta + 2 - iv)
      zcom(iv) =-zcom(nzeta + 2 - iv)

input.cma reaches guess_axis in both symmetry modes. With LASYM = F the reguess succeeds and the run goes on to ns = 51. With LASYM = T and every asymmetric coefficient left at zero, so the same equilibrium, it fails:

  NS =   25 NO. FOURIER MODES =   59 FTOLV =  1.000E-06 NITER =   1000
  INITIAL JACOBIAN CHANGED SIGN!
  TRYING TO IMPROVE INITIAL MAGNETIC AXIS GUESS

 INITIAL JACOBIAN CHANGED SIGN (IMPROVE INITIAL GUESS)

 The jacobian was non-definite!

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions