Initialize the coarse retry after a bad initial Jacobian

src/vmec.f90

@@ -42,6 +42,8 @@ subroutine vmec
   INTEGER :: ns_old=0
   INTEGER :: igrid
   INTEGER :: jacob_off
+  REAL(rprec), DIMENSION(0:ntord) :: initial_raxis_cc, initial_raxis_cs
+  REAL(rprec), DIMENSION(0:ntord) :: initial_zaxis_cc, initial_zaxis_cs
 
   integer :: i, js
 
@@ -99,6 +101,11 @@ subroutine vmec
      ! COMPUTE NS-INVARIANT ARRAYS
      CALL fixaray
 
+     initial_raxis_cc = raxis_cc
+     initial_raxis_cs = raxis_cs
+     initial_zaxis_cc = zaxis_cc
+     initial_zaxis_cs = zaxis_cs
+
      ns_old = 0
      delt0r = delt
 
@@ -199,11 +206,19 @@ subroutine vmec
         ! if did not converge only because jacobian was bad
         ! and the intermediate ns=3 run was not performed yet (jacob_off is still == 0),
         ! retry the whole thing again
-        IF (ier_flag.ne.bad_jacobian_flag) THEN
+        IF (ier_flag.ne.bad_jacobian_flag .or. jacob_off.ne.0 .or. &
+            igrid.ne.1 .or. ns_array(1).le.3) THEN
            exit ! jacob_off loop
            ! otherwise, retry with initial ns=3 to fix bad jacobian
         END IF
 
+        ! Initialize a cold three-surface mesh from the input axis.
+        ns_old = 0
+        raxis_cc = initial_raxis_cc
+        raxis_cs = initial_raxis_cs
+        zaxis_cc = initial_zaxis_cc
+        zaxis_cs = initial_zaxis_cs
+
      ! if ier_flag .eq. bad_jacobian_flag, repeat once again with ns=3 before
      end do ! jacob_off = 0, 1
 

test/unit/CMakeLists.txt

@@ -1,6 +1,6 @@
-# Unit tests for individual numerical routines.
+# Numerical unit tests and a solver-level initialization regression.
 #
-# Standalone targets that compile only the sources under test, so they do not
+# The spline target compiles only the sources under test, so it does not
 # depend on the full vmec library or its NetCDF/FFTW/LAPACK link.
 
 add_executable (test_spline_akima
@@ -10,3 +10,17 @@ add_executable (test_spline_akima
 set_target_properties (test_spline_akima PROPERTIES
   Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/modules)
 add_test (NAME spline_akima_reflection_symmetry COMMAND test_spline_akima)
+
+add_executable(check_jacobian_recovery
+  check_jacobian_recovery.f90
+  ${CMAKE_SOURCE_DIR}/src/data/stel_kinds.f90)
+set_target_properties(check_jacobian_recovery PROPERTIES
+  Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/recovery_modules)
+target_link_libraries(check_jacobian_recovery PRIVATE NetCDF::NetCDF)
+add_test(NAME initial_jacobian_recovery
+  COMMAND ${CMAKE_COMMAND}
+    -DVMEC=$<TARGET_FILE:xvmec>
+    -DCHECKER=$<TARGET_FILE:check_jacobian_recovery>
+    -DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/input.initial_jacobian_recovery
+    -DWORK_DIR=${CMAKE_CURRENT_BINARY_DIR}/jacobian_recovery
+    -P ${CMAKE_CURRENT_SOURCE_DIR}/run_jacobian_recovery.cmake)

test/unit/check_jacobian_recovery.f90

@@ -0,0 +1,38 @@
+program check_jacobian_recovery
+  use netcdf
+  use stel_kinds, only: rprec
+  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
+  implicit none
+
+  integer :: ncid, varid, ns, ier, i
+  real(rprec) :: residuals(3)
+  character(len=4), parameter :: names(3) = [character(len=4) :: 'fsqr', 'fsqz', 'fsql']
+  character(len=512) :: filename
+
+  call get_command_argument(1, filename)
+  call checked(nf90_open(trim(filename), nf90_nowrite, ncid))
+  call checked(nf90_inq_varid(ncid, 'ns', varid))
+  call checked(nf90_get_var(ncid, varid, ns))
+  call checked(nf90_inq_varid(ncid, 'ier_flag', varid))
+  call checked(nf90_get_var(ncid, varid, ier))
+  do i = 1, 3
+    call checked(nf90_inq_varid(ncid, names(i), varid))
+    call checked(nf90_get_var(ncid, varid, residuals(i)))
+  end do
+  call checked(nf90_close(ncid))
+
+  if (ns /= 31 .or. ier /= 0) error stop 'did not reach the requested equilibrium'
+  if (.not. all(ieee_is_finite(residuals))) error stop 'non-finite force residual'
+  if (any(residuals > 1.0e-9_rprec)) error stop 'requested force tolerance was not met'
+
+contains
+
+  subroutine checked(status)
+    integer, intent(in) :: status
+    if (status /= nf90_noerr) then
+      print *, trim(nf90_strerror(status))
+      error stop 'cannot read equilibrium output'
+    end if
+  end subroutine checked
+
+end program check_jacobian_recovery

test/unit/initial_jacobian_recovery.md

@@ -0,0 +1,10 @@
+`input.initial_jacobian_recovery` imposes the outermost boundary of
+[QUASR configuration 14802](https://quasr.flatironinstitute.org/simsopt_serials/0014/serial0014802.json)
+at `mpol=ntor=6`, with zero pressure and prescribed zero toroidal current.
+The toroidal flux comes from the prescribed coil vector potential. Its boundary
+coefficients were converted with SIMSOPT and are not changed during recovery.
+
+The cold `[8,16,31]` sequence has a bad initial Jacobian after the axis search.
+An explicit `[3,8,16,31]` sequence converges, using `ftol=1e-4` only for the
+three-surface bootstrap and `1e-9` for the requested grids. The regression test
+requires the automatic recovery to reach the requested `ns=31` and tolerance.

test/unit/input.initial_jacobian_recovery

@@ -0,0 +1,175 @@
+&INDATA
+ lasym = .false.
+ nfp = 2
+ mpol = 6
+ ntor = 6
+ ntheta = 0
+ nzeta = 24
+ ns_array = 8, 16, 31
+ ftol_array = 1e-09, 1e-09, 1e-09
+ niter_array = 4000, 4000, 4000
+ phiedge = 0.17677284194671744
+ ncurr = 1
+ pmass_type = 'power_series'
+ am = 1.0, -1.0
+ pres_scale = 0.0
+ gamma = 0.0
+ spres_ped = 1.0
+ piota_type = 'power_series'
+ pcurr_type = 'power_series'
+ ac = 1.0, -1.0
+ curtor = 0.0
+ bloat = 1.0
+ lfreeb = .false.
+ mgrid_file = 'NONE'
+ nvacskip = 6
+ nstep = 200
+ aphi = 1.0
+ delt = 0.25
+ tcon0 = 1.0
+ raxis_cc = 0.9929268048162759, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+ zaxis_cs = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+ rbc(0,0) = 0.9929268048162759
+ rbc(1,0) = 0.15683336219439925
+ rbc(2,0) = 0.008518638515756934
+ rbc(3,0) = -0.0006540562935552495
+ rbc(4,0) = 7.375680186709275e-05
+ rbc(5,0) = 0.00029873994273236403
+ rbc(6,0) = 0.00014028054389655275
+ rbc(-6,1) = -7.816701625651354e-05
+ rbc(-5,1) = -3.486321790835295e-05
+ rbc(-4,1) = 0.00012620428698713207
+ rbc(-3,1) = -0.0004669922867383307
+ rbc(-2,1) = -0.009120850549362116
+ rbc(-1,1) = -0.054328197386995895
+ rbc(0,1) = -0.2835163396789723
+ rbc(1,1) = 0.1587135742967791
+ rbc(2,1) = 0.03551477617973609
+ rbc(3,1) = 0.004441191341754036
+ rbc(4,1) = 0.0003845034555506475
+ rbc(5,1) = 5.8299922125630376e-05
+ rbc(6,1) = 0.00012819278728035322
+ rbc(-6,2) = -1.2970119548189213e-05
+ rbc(-5,2) = -7.476499602708074e-05
+ rbc(-4,2) = -0.0001585781696243108
+ rbc(-3,2) = -0.00010210370651865964
+ rbc(-2,2) = 0.0014840022354688748
+ rbc(-1,2) = 0.0028683637734258334
+ rbc(0,2) = -0.013716269081447992
+ rbc(1,2) = 0.04779887964954569
+ rbc(2,2) = 0.03402446603211992
+ rbc(3,2) = 0.008768080620779495
+ rbc(4,2) = 0.0012503442717795498
+ rbc(5,2) = -1.9765179563693993e-05
+ rbc(6,2) = -7.223333507366949e-05
+ rbc(-6,3) = 9.837530027517543e-06
+ rbc(-5,3) = 1.2932019954255361e-05
+ rbc(-4,3) = -4.169904825704345e-05
+ rbc(-3,3) = -0.00035981922821330225
+ rbc(-2,3) = -0.001118593982795612
+ rbc(-1,3) = 0.0008855817417740778
+ rbc(0,3) = -0.000650942355837123
+ rbc(1,3) = 0.006725470426319122
+ rbc(2,3) = 0.016330639060982716
+ rbc(3,3) = 0.008592600755727932
+ rbc(4,3) = 0.0023189155652774563
+ rbc(5,3) = 0.00035720662074978
+ rbc(6,3) = -6.85835987172996e-05
+ rbc(-6,4) = -1.0229376954031479e-05
+ rbc(-5,4) = -1.450057112686764e-05
+ rbc(-4,4) = -1.6519906000237664e-05
+ rbc(-3,4) = -4.7203955436374055e-05
+ rbc(-2,4) = -0.00027934041253135466
+ rbc(-1,4) = 0.000253627523048619
+ rbc(0,4) = -0.0002675889437167926
+ rbc(1,4) = 0.0020832051008993304
+ rbc(2,4) = 0.0054818794515777486
+ rbc(3,4) = 0.005288710591377669
+ rbc(4,4) = 0.0024153835165718167
+ rbc(5,4) = 0.0007108759652516153
+ rbc(6,4) = 0.00011517199173844234
+ rbc(-6,5) = 3.6248771095641698e-06
+ rbc(-5,5) = 9.503994063217063e-06
+ rbc(-4,5) = 1.5916300686587047e-05
+ rbc(-3,5) = 3.036231525323081e-05
+ rbc(-2,5) = 4.6678036831549606e-05
+ rbc(-1,5) = 0.000373147088586584
+ rbc(0,5) = 0.0008852005721704652
+ rbc(1,5) = 0.0003387808897669372
+ rbc(2,5) = 0.001193348135974064
+ rbc(3,5) = 0.0022281821286734486
+ rbc(4,5) = 0.0017414695876251094
+ rbc(5,5) = 0.0007810763137289322
+ rbc(6,5) = 0.0002569673010588241
+ zbs(1,0) = 0.1305846106532369
+ zbs(2,0) = 0.00808804122597515
+ zbs(3,0) = -0.0005609874006237156
+ zbs(4,0) = -0.0004638755241490486
+ zbs(5,0) = -5.792340604349128e-05
+ zbs(6,0) = 5.723499427861041e-05
+ zbs(-6,1) = -1.1925002821131885e-05
+ zbs(-5,1) = -1.3721596734418438e-05
+ zbs(-4,1) = 0.00014424033909440776
+ zbs(-3,1) = 0.0012245136833539064
+ zbs(-2,1) = 0.00672527120922041
+ zbs(-1,1) = 0.048372018802395216
+ zbs(0,1) = 0.31334019322069073
+ zbs(1,1) = 0.1302309561875805
+ zbs(2,1) = 0.0317897924579535
+ zbs(3,1) = 0.0035891116134236845
+ zbs(4,1) = 0.000488650219018151
+ zbs(5,1) = -1.6760775894436837e-05
+ zbs(6,1) = -1.3831257409952964e-05
+ zbs(-6,2) = -1.814716612259062e-05
+ zbs(-5,2) = -6.55651725966855e-05
+ zbs(-4,2) = -9.247032481445009e-05
+ zbs(-3,2) = -0.0003495262254473919
+ zbs(-2,2) = -0.0016144345817882219
+ zbs(-1,2) = -0.004072385342281151
+ zbs(0,2) = 0.0009446954755320977
+ zbs(1,2) = 0.03765698634784408
+ zbs(2,2) = 0.031539227315748486
+ zbs(3,2) = 0.008055208750449275
+ zbs(4,2) = 0.0012486811294420482
+ zbs(5,2) = 0.00016950400247146068
+ zbs(6,2) = -3.839046405268226e-05
+ zbs(-6,3) = 2.3384711634698186e-05
+ zbs(-5,3) = 3.157912257933612e-05
+ zbs(-4,3) = 3.621718480354919e-06
+ zbs(-3,3) = -1.5892234933467628e-05
+ zbs(-2,3) = -0.0003307950411325112
+ zbs(-1,3) = 0.0007680049021077882
+ zbs(0,3) = -0.003144595013522224
+ zbs(1,3) = 0.004266680507001436
+ zbs(2,3) = 0.014861102355524496
+ zbs(3,3) = 0.008422018761421925
+ zbs(4,3) = 0.00225604105407406
+ zbs(5,3) = 0.0004262121224184639
+ zbs(6,3) = 3.658292794565433e-05
+ zbs(-6,4) = 9.39604698317768e-06
+ zbs(-5,4) = 1.859949791748088e-05
+ zbs(-4,4) = 2.072215627624109e-05
+ zbs(-3,4) = 8.278195748988871e-06
+ zbs(-2,4) = -0.00019510476317084475
+ zbs(-1,4) = 0.0005750361234318032
+ zbs(0,4) = 0.0008582765616717533
+ zbs(1,4) = 0.002691336474273391
+ zbs(2,4) = 0.004778968386080365
+ zbs(3,4) = 0.005132454324468034
+ zbs(4,4) = 0.0024565184347139345
+ zbs(5,4) = 0.0007038908044132183
+ zbs(6,4) = 0.0001433335639518173
+ zbs(-6,5) = -2.1193847038565956e-06
+ zbs(-5,5) = -1.6572873451586166e-06
+ zbs(-4,5) = 1.641127551435855e-05
+ zbs(-3,5) = 4.8137690542497766e-05
+ zbs(-2,5) = -5.945535837675628e-07
+ zbs(-1,5) = 3.1341902479503886e-06
+ zbs(0,5) = 2.2900941016509426e-05
+ zbs(1,5) = -0.0009793432199695416
+ zbs(2,5) = 0.0009375747360851498
+ zbs(3,5) = 0.0021209886671250473
+ zbs(4,5) = 0.0017332101629774824
+ zbs(5,5) = 0.0007834491521010137
+ zbs(6,5) = 0.0002445886962532652
+/

test/unit/run_jacobian_recovery.cmake

@@ -0,0 +1,25 @@
+file(MAKE_DIRECTORY "${WORK_DIR}")
+configure_file("${INPUT}" "${WORK_DIR}/input.recovery" COPYONLY)
+file(REMOVE "${WORK_DIR}/wout_recovery.nc")
+execute_process(
+  COMMAND "${VMEC}" input.recovery
+  WORKING_DIRECTORY "${WORK_DIR}"
+  RESULT_VARIABLE result
+  OUTPUT_FILE "${WORK_DIR}/solver.log"
+  ERROR_FILE "${WORK_DIR}/solver.err"
+  TIMEOUT 120)
+if(NOT result EQUAL 0)
+  message(FATAL_ERROR "VMEC failed: ${result}")
+endif()
+file(READ "${WORK_DIR}/solver.log" log)
+string(REGEX MATCHALL "NS = +3 +NO[.] FOURIER" coarse_stages "${log}")
+list(LENGTH coarse_stages coarse_count)
+if(NOT coarse_count EQUAL 1)
+  message(FATAL_ERROR "Expected one three-surface recovery stage, got ${coarse_count}")
+endif()
+execute_process(
+  COMMAND "${CHECKER}" "${WORK_DIR}/wout_recovery.nc"
+  RESULT_VARIABLE result)
+if(NOT result EQUAL 0)
+  message(FATAL_ERROR "The recovered equilibrium failed validation")
+endif()