Synchronize periodic multispecies refinement

NEO-2-QL/CMakeSources.in

@@ -6,6 +6,7 @@ set(NEO2_QL_SRC_FILES
   helical_response_mod.f90
   qflux_profile_mod.f90
   phi_crossing_alignment_mod.f90
+  multispecies_iteration_sync_mod.f90
   propagator.f90
   flint.f90
   ripple_solver_axi_test.f90

NEO-2-QL/multispecies_iteration_sync_mod.f90

@@ -0,0 +1,35 @@
+MODULE multispecies_iteration_sync_mod
+    USE mpiprovider_module, ONLY : mpro
+    IMPLICIT NONE
+    PRIVATE
+    PUBLIC :: gather_iteration_convergence
+
+CONTAINS
+
+    ! Every integral_part call contains a world allgather.  Species may have
+    ! different local convergence rates, so all ranks must continue the outer
+    ! iteration until every species has converged; otherwise a fast species can
+    ! enter the next propagator collective while a slow species is still here.
+    SUBROUTINE gather_iteration_convergence(local_change,local_threshold, &
+            ispec,changes,thresholds,converged)
+        DOUBLE PRECISION, INTENT(IN) :: local_change,local_threshold
+        INTEGER, INTENT(IN) :: ispec
+        DOUBLE PRECISION, DIMENSION(0:), INTENT(INOUT) :: changes,thresholds
+        LOGICAL, INTENT(OUT) :: converged
+
+        IF (SIZE(changes) .NE. mpro%getNumProcs() .OR. &
+            SIZE(thresholds) .NE. mpro%getNumProcs()) THEN
+            ERROR STOP 'multispecies convergence arrays do not match MPI size'
+        END IF
+        IF (ispec .LT. 0 .OR. ispec .GE. mpro%getNumProcs()) THEN
+            ERROR STOP 'multispecies convergence species index is out of range'
+        END IF
+
+        changes(ispec)=local_change
+        thresholds(ispec)=local_threshold
+        CALL mpro%allgather_inplace(changes)
+        CALL mpro%allgather_inplace(thresholds)
+        converged=ALL(changes .LT. thresholds)
+    END SUBROUTINE gather_iteration_convergence
+
+END MODULE multispecies_iteration_sync_mod

NEO-2-QL/propagator.f90

@@ -1672,6 +1672,12 @@ END SUBROUTINE ripple_solver_ArnoldiO2
     ENDIF
     !! End Modification by Andreas F. Martitsch (27.07.2015)
 
+    ! A recoverable crossing/refinement request (ierr=3) returns before the
+    ! ripple solver allocates the propagator matrices.  The caller handles that
+    ! status by refining phi and retrying; do not dereference the deliberately
+    ! absent matrices on this error path.
+    IF (ierr .NE. 0) RETURN
+
     prop_c%p%npass_l = SIZE(prop_c%p%amat_p_m,1)/(prop_c%p%nvelocity+1)
     prop_c%p%npass_r = SIZE(prop_c%p%amat_m_p,1)/(prop_c%p%nvelocity+1)
 

NEO-2-QL/qflux_profile_mod.f90

@@ -109,23 +109,27 @@ subroutine qflux_point_components_from_flux_vector(flux_row, source_col, &
         end do
     end subroutine qflux_point_components_from_flux_vector
 
-    subroutine record_qflux_interface_traces(tag, phi, eta_grid, block_base, &
+    subroutine record_qflux_interface_traces(tag, phi, bhat, eta_grid, block_base, &
             block_npassing, lag, step_plus, step_minus, flux_row, source_rhs, &
             source_solution, ierr)
         integer, intent(in) :: tag, block_base(:), block_npassing(:), lag
-        real(dp), intent(in) :: phi(:), eta_grid(0:), step_plus(:), step_minus(:)
+        real(dp), intent(in) :: phi(:), bhat(:), eta_grid(0:)
+        real(dp), intent(in) :: step_plus(:), step_minus(:)
         real(dp), intent(in) :: flux_row(:), source_rhs(:, :)
         real(dp), intent(in) :: source_solution(:, :)
         integer, intent(out) :: ierr
         character(len=5) :: side
         integer :: band, base, column, endpoint, force, iunit, m, nband, point
         integer :: status
-        real(dp) :: constant_coordinate(3), source_solution_orthogonal
+        real(dp), allocatable :: null_state(:)
+        real(dp) :: constant_coordinate(3), null_measure
+        real(dp) :: source_solution_orthogonal
 
         ierr = 0
         if (.not. interface_trace_initialized) call initialize_interface_trace(ierr)
         if (ierr /= 0 .or. .not. allocated(interface_trace_filename)) return
         if (size(phi) /= size(block_base) &
+            .or. size(bhat) /= size(block_base) &
             .or. size(block_npassing) /= size(block_base) &
             .or. size(step_plus) /= size(block_base) &
             .or. size(step_minus) /= size(block_base) &
@@ -139,12 +143,14 @@ subroutine record_qflux_interface_traces(tag, phi, eta_grid, block_base, &
         if (size(block_base) < 2 &
             .or. maxval(block_npassing) >= size(eta_grid) &
             .or. any(block_npassing < 0) &
+            .or. any(bhat <= 0.0_dp) &
             .or. any(step_plus == 0.0_dp) &
             .or. any(step_minus == 0.0_dp)) then
             ierr = 1
             return
         end if
         if (.not. all(ieee_is_finite(phi)) &
+            .or. .not. all(ieee_is_finite(bhat)) &
             .or. .not. all(ieee_is_finite(eta_grid)) &
             .or. .not. all(ieee_is_finite(step_plus)) &
             .or. .not. all(ieee_is_finite(step_minus)) &
@@ -171,17 +177,27 @@ subroutine record_qflux_interface_traces(tag, phi, eta_grid, block_base, &
             end if
             nband = block_npassing(point) + 1
             base = block_base(point)
-            ! The Lorentz collision operator's discrete right-null state is
-            ! one constant shared by the co- and counter-passing lag=0 rows.
-            ! Export the Euclidean-orthogonal component as a diagnostic; the
-            ! raw solution and all solver inputs remain unchanged.
+            ! In the finite-volume representation a constant distribution is
+            ! stored as the pitch-bin width vector, not as a vector of ones.
+            ! Remove its mass coordinate, shared by the co- and counter-
+            ! passing lag=0 rows.  This is the same left/right null pairing
+            ! used by the Lorentz solvability projection in join_ends.
+            allocate(null_state(nband))
+            if (nband > 1) null_state(1:nband - 1) = &
+                eta_grid(1:nband - 1) - eta_grid(0:nband - 2)
+            null_state(nband) = 1.0_dp/bhat(point) - eta_grid(nband - 1)
+            if (any(null_state <= 0.0_dp)) then
+                status = 1
+                exit
+            end if
+            null_measure = 2.0_dp*sum(null_state)
             do force = 1, 3
                 constant_coordinate(force) = sum( &
                     source_solution(base + 1:base + nband, force))
                 constant_coordinate(force) = constant_coordinate(force) + sum( &
                     source_solution(base + nband + 1:base + 2*nband, force))
             end do
-            constant_coordinate = constant_coordinate/real(2*nband, dp)
+            constant_coordinate = constant_coordinate/null_measure
             do m = 0, lag
                 base = block_base(point) + 2*m*nband
                 do band = 1, nband
@@ -191,7 +207,7 @@ subroutine record_qflux_interface_traces(tag, phi, eta_grid, block_base, &
                             source_solution(column, force)
                         if (m == 0) source_solution_orthogonal = &
                             source_solution_orthogonal - &
-                            constant_coordinate(force)
+                            null_state(band)*constant_coordinate(force)
                         call write_interface_trace(iunit, tag, side, 'p', m, &
                             force, band - 1, eta_grid(band - 1), phi(point), &
                             flux_row(column)/step_plus(point), &
@@ -207,7 +223,7 @@ subroutine record_qflux_interface_traces(tag, phi, eta_grid, block_base, &
                             source_solution(column, force)
                         if (m == 0) source_solution_orthogonal = &
                             source_solution_orthogonal - &
-                            constant_coordinate(force)
+                            null_state(band)*constant_coordinate(force)
                         call write_interface_trace(iunit, tag, side, 'm', m, &
                             force, band - 1, eta_grid(band - 1), phi(point), &
                             flux_row(column)/step_minus(point), &
@@ -220,6 +236,7 @@ subroutine record_qflux_interface_traces(tag, phi, eta_grid, block_base, &
                 end do
                 if (status /= 0) exit
             end do
+            deallocate(null_state)
             if (status /= 0) exit
         end do
         close(iunit, iostat=ierr)

NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90

@@ -86,6 +86,8 @@ subroutine ripple_solver_ArnoldiO2( &
     ! MPI SUPPORT for multi-species part
     ! (run with, e.g.,  mpiexec -np 3 ./neo2.x)
     use mpiprovider_module, only: mpro
+    use mpi, only: MPI_ALLREDUCE, MPI_COMM_WORLD, MPI_INTEGER, MPI_MAX, &
+                   MPI_SUCCESS
     ! Load x1mm and x2mm (=energy dependence of drift frequencies)
     ! from collision operator module. This step allows for
     ! support of different basis functions and replaces routine "lagxmm".
@@ -277,6 +279,7 @@ subroutine ripple_solver_ArnoldiO2( &
   !! Modification by Andreas F. Martitsch (28.07.2015)
     !  multi-species part
     integer :: ispec, ispecp, ispecpp ! species indices
+    integer :: refine_local, refine_global, mpi_ierr
     integer :: drive_spec
     complex(dp), dimension(:, :, :), allocatable :: source_vector_all
     real(dp), dimension(:, :, :, :), allocatable :: qflux_allspec
@@ -825,7 +828,14 @@ subroutine ripple_solver_ArnoldiO2( &
     geodcu_mfl = arr_comp(:, 1)
     deallocate (arr_real, arr_comp)
 
-    if (maxval(phi_divide) > 1) then
+    refine_local = merge(1, 0, maxval(phi_divide) > 1)
+    call MPI_ALLREDUCE(refine_local, refine_global, 1, MPI_INTEGER, MPI_MAX, &
+                       MPI_COMM_WORLD, mpi_ierr)
+    if (mpi_ierr /= MPI_SUCCESS) then
+        ierr = 1
+        return
+    end if
+    if (refine_global /= 0) then
         ierr = 3
         write (*, *) 'ERROR: crossing geometry requires phi refinement.'
         return

NEO-2-QL/ripple_solver_axi_test.f90

@@ -86,6 +86,9 @@ SUBROUTINE ripple_solver(                                 &
                             nz_regper,irow_regper,icol_regper,amat_regper,     &
                             deallocate_ntv_eqmat
   use mpiprovider_module, only : mpro
+  USE multispecies_iteration_sync_mod, ONLY : gather_iteration_convergence
+  USE mpi, ONLY : MPI_ALLREDUCE, MPI_COMM_WORLD, MPI_INTEGER, MPI_MAX, &
+       MPI_SUCCESS
   USE collop
   USE helical_source_mod, ONLY : add_helical_source, &
        apply_reconstructed_incoming_rows
@@ -276,6 +279,7 @@ SUBROUTINE ripple_solver(                                 &
   !! Modification by Andreas F. Martitsch (28.07.2015)
   !  multi-species part
   INTEGER :: ispec, ispecp, ispecpp ! species indices
+  INTEGER :: refine_local,refine_global,mpi_ierr
   INTEGER :: crossing_ncomp, crossing_nreal
   INTEGER :: drive_spec
   INTEGER :: isw_regper, ipart1
@@ -286,6 +290,7 @@ SUBROUTINE ripple_solver(                                 &
   DOUBLE PRECISION, DIMENSION(0:num_spec-1) :: break_cond1
   DOUBLE PRECISION, DIMENSION(0:num_spec-1) :: break_cond2
   INTEGER :: prop_fileformat=0
+  LOGICAL :: iteration_converged
   !! End Modification by Andreas F. Martitsch (28.07.2015)
 
   !***************************
@@ -800,7 +805,14 @@ SUBROUTINE ripple_solver(                                 &
   geodcu_mfl=real(crossing_complex(:,1),kind=dp)
   deallocate(crossing_real,crossing_complex,unused_dbcovar_mfl)
 
-  if(maxval(phi_divide).gt.1) then
+  refine_local=MERGE(1,0,maxval(phi_divide).gt.1)
+  CALL MPI_ALLREDUCE(refine_local,refine_global,1,MPI_INTEGER,MPI_MAX, &
+       MPI_COMM_WORLD,mpi_ierr)
+  IF (mpi_ierr .NE. MPI_SUCCESS) THEN
+    ierr=1
+    RETURN
+  END IF
+  if(refine_global.ne.0) then
     ierr=3
     write(*,*) 'ERROR: crossing geometry requires phi refinement.'
     return
@@ -2947,7 +2959,8 @@ SUBROUTINE ripple_solver(                                 &
 IF (isw_hel_drive.NE.0 .AND. num_spec.EQ.1 .AND. mpro%getrank().EQ.0) THEN
    IF (iplot.EQ.1) THEN
       CALL record_qflux_interface_traces(fieldpropagator%tag, &
-           phi_mfl(ibeg:iend),eta,ind_start(ibeg:iend),npl(ibeg:iend),lag, &
+           phi_mfl(ibeg:iend),bhat_mfl(ibeg:iend),eta, &
+           ind_start(ibeg:iend),npl(ibeg:iend),lag, &
            cp_step_p,cp_step_m,flux_vector(2,:),source_vector, &
            source_vector_all(:,:,0),cp_trace_status)
       IF (cp_trace_status.NE.0) &
@@ -3125,10 +3138,12 @@ SUBROUTINE ripple_solver(                                 &
 
           bvec_sp=bvec_lor+bvec_iter
 
-          if(sum(abs(bvec_sp-bvec_prev)) .lt.                                 &
-             sum(abs(bvec_prev))*epserr_iter) then
-            exit
-          endif
+          ! integral_part uses a world allgather; keep every species on the
+          ! same iteration count before the next propagator collective.
+          CALL gather_iteration_convergence( &
+               SUM(ABS(bvec_sp-bvec_prev)),SUM(ABS(bvec_prev))*epserr_iter, &
+               ispec,break_cond1,break_cond2,iteration_converged)
+          IF(iteration_converged) EXIT
 
         enddo
 
@@ -3181,10 +3196,10 @@ SUBROUTINE ripple_solver(                                 &
 
           bvec_sp=bvec_lor+bvec_iter
 
-          if(sum(abs(bvec_sp-bvec_prev)) .lt.                                 &
-             sum(abs(bvec_prev))*epserr_iter) then
-             exit
-          endif
+          CALL gather_iteration_convergence( &
+               SUM(ABS(bvec_sp-bvec_prev)),SUM(ABS(bvec_prev))*epserr_iter, &
+               ispec,break_cond1,break_cond2,iteration_converged)
+          IF(iteration_converged) EXIT
 
         enddo
 

TEST/CMakeLists.txt

@@ -234,6 +234,23 @@ set_tests_properties(ripple_solver_route_test PROPERTIES
     TIMEOUT 30
 )
 
+add_executable(test_multispecies_iteration_sync
+    test_multispecies_iteration_sync.f90)
+target_link_libraries(test_multispecies_iteration_sync neo2_ql)
+
+add_test(NAME multispecies_iteration_sync_test
+         COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2
+                 ${MPIEXEC_PREFLAGS}
+                 $<TARGET_FILE:test_multispecies_iteration_sync>
+                 ${MPIEXEC_POSTFLAGS})
+
+set_tests_properties(multispecies_iteration_sync_test PROPERTIES
+    PASS_REGULAR_EXPRESSION "All tests passed!"
+    FAIL_REGULAR_EXPRESSION "FAIL"
+    PROCESSORS 2
+    TIMEOUT 30
+)
+
 add_executable(test_helical_source test_helical_source.f90)
 target_link_libraries(test_helical_source neo2_ql)
 

TEST/test_multispecies_iteration_sync.f90

@@ -0,0 +1,43 @@
+PROGRAM test_multispecies_iteration_sync
+    USE mpiprovider_module, ONLY : mpro
+    USE multispecies_iteration_sync_mod, ONLY : gather_iteration_convergence
+    IMPLICIT NONE
+
+    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: changes,thresholds
+    DOUBLE PRECISION :: local_change
+    INTEGER :: ispec,iter,exit_iteration
+    LOGICAL :: converged
+
+    CALL mpro%init()
+    IF (mpro%getNumProcs() .NE. 2) &
+        ERROR STOP 'FAIL: test requires exactly two MPI ranks'
+    ispec=mpro%getRank()
+    ALLOCATE(changes(0:mpro%getNumProcs()-1))
+    ALLOCATE(thresholds(0:mpro%getNumProcs()-1))
+    changes=0.d0
+    thresholds=0.d0
+    exit_iteration=0
+
+    DO iter=1,5
+        IF (iter .GE. ispec+2) THEN
+            local_change=0.25d0
+        ELSE
+            local_change=1.d0
+        END IF
+        CALL gather_iteration_convergence(local_change,0.5d0,ispec, &
+            changes,thresholds,converged)
+        IF (converged) THEN
+            exit_iteration=iter
+            EXIT
+        END IF
+    END DO
+
+    IF (exit_iteration .NE. 3) &
+        ERROR STOP 'FAIL: ranks did not wait for the slow species'
+    IF (ANY(changes .GT. thresholds)) &
+        ERROR STOP 'FAIL: gathered convergence state is inconsistent'
+    IF (mpro%isMaster()) PRINT *, 'All tests passed!'
+
+    DEALLOCATE(changes,thresholds)
+    CALL mpro%deinit(.FALSE.)
+END PROGRAM test_multispecies_iteration_sync