From b653480c3382721624e8ad3d1e780a1b17832434 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 21 Jul 2026 10:35:59 +0200 Subject: [PATCH 1/5] Project the finite-volume Lorentz null coordinate --- NEO-2-QL/qflux_profile_mod.f90 | 37 +++++++++++++++++------- NEO-2-QL/ripple_solver_axi_test.f90 | 3 +- TEST/test_qflux_interface_diagnostic.f90 | 24 +++++++++++---- 3 files changed, 48 insertions(+), 16 deletions(-) diff --git a/NEO-2-QL/qflux_profile_mod.f90 b/NEO-2-QL/qflux_profile_mod.f90 index de3b09c5..834bf837 100644 --- a/NEO-2-QL/qflux_profile_mod.f90 +++ b/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) diff --git a/NEO-2-QL/ripple_solver_axi_test.f90 b/NEO-2-QL/ripple_solver_axi_test.f90 index e3c34057..2bdad1c3 100644 --- a/NEO-2-QL/ripple_solver_axi_test.f90 +++ b/NEO-2-QL/ripple_solver_axi_test.f90 @@ -2947,7 +2947,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) & diff --git a/TEST/test_qflux_interface_diagnostic.f90 b/TEST/test_qflux_interface_diagnostic.f90 index 826d6f9e..be52ddd8 100644 --- a/TEST/test_qflux_interface_diagnostic.f90 +++ b/TEST/test_qflux_interface_diagnostic.f90 @@ -6,16 +6,19 @@ program test_qflux_interface_diagnostic character(len=512) :: filename, header character(len=5) :: side character(len=1) :: direction - integer :: band, force, i, ierr, iunit, laguerre, sequence, status, tag + integer :: band, endpoint, force, i, ierr, iunit, laguerre, sequence + integer :: status, tag integer :: block_base(2), block_npassing(2), rows - real(real64) :: eta, eta_grid(0:1), flux_kernel, flux_row(8), phi + real(real64) :: bhat(2), eta, eta_grid(0:1), flux_kernel, flux_row(8), phi real(real64) :: phi_grid(2), solution(8, 3), source_rhs, rhs(8, 3) + real(real64) :: projected_sum(2, 3) real(real64) :: projected, step_minus(2), step_plus(2), value logical :: found_counter_record block_base = [0, 4] block_npassing = 1 - eta_grid = [0.25_real64, 0.75_real64] + bhat = [1.0_real64, 0.8_real64] + eta_grid = [0.0_real64, 0.75_real64] phi_grid = [1.0_real64, 2.0_real64] step_plus = [2.0_real64, 4.0_real64] step_minus = [5.0_real64, 10.0_real64] @@ -24,7 +27,7 @@ program test_qflux_interface_diagnostic shape(rhs)) solution = rhs + 1000.0_real64 - call record_qflux_interface_traces(12, phi_grid, eta_grid, block_base, & + call record_qflux_interface_traces(12, phi_grid, bhat, eta_grid, block_base, & block_npassing, 0, step_plus, step_minus, flux_row, rhs, solution, ierr) if (ierr /= 0) error stop 'FAIL: interface trace record returned an error' call get_environment_variable('NEO2_INTERFACE_TRACE_FILE', filename) @@ -40,6 +43,7 @@ program test_qflux_interface_diagnostic rows = 0 found_counter_record = .false. + projected_sum = 0.0_real64 do read(iunit, *, iostat=status) sequence, tag, side, direction, laguerre, & force, band, eta, phi, flux_kernel, source_rhs, value, projected @@ -48,6 +52,14 @@ program test_qflux_interface_diagnostic rows = rows + 1 if (sequence /= rows .or. tag /= 12 .or. laguerre /= 0) & error stop 'FAIL: interface trace identity is wrong' + if (trim(side) == 'left') then + endpoint = 1 + else if (trim(side) == 'right') then + endpoint = 2 + else + error stop 'FAIL: interface trace side is wrong' + end if + projected_sum(endpoint, force) = projected_sum(endpoint, force) + projected if (trim(side) == 'right' .and. direction == 'm' .and. band == 1 & .and. force == 3) then found_counter_record = .true. @@ -57,13 +69,15 @@ program test_qflux_interface_diagnostic .or. abs(source_rhs - rhs(7, 3)) > 1.0e-15_real64 & .or. abs(value - solution(7, 3)) > 1.0e-15_real64 & .or. abs(projected - (solution(7, 3) & - - sum(solution(5:8, 3))/4.0_real64)) & + - 0.5_real64*sum(solution(5:8, 3))/2.5_real64)) & > 1.0e-15_real64) & error stop 'FAIL: counter-passing trace mapping is wrong' end if end do close(iunit) if (rows /= 24) error stop 'FAIL: interface trace row count is wrong' + if (maxval(abs(projected_sum)) > 1.0e-10_real64) & + error stop 'FAIL: projected trace retains the finite-volume null mass' if (.not. found_counter_record) & error stop 'FAIL: expected counter-passing trace record is absent' From d8e765b0ed6b95f78e85ab97c7f0da74c6f0acd6 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 21 Jul 2026 11:46:58 +0200 Subject: [PATCH 2/5] Honor recoverable ripple-solver refinement status --- NEO-2-QL/propagator.f90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEO-2-QL/propagator.f90 b/NEO-2-QL/propagator.f90 index 76f8f033..3956400f 100644 --- a/NEO-2-QL/propagator.f90 +++ b/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) From ae418a77a7df8d12e65dbe475371587d3ad7f0f7 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 21 Jul 2026 12:00:59 +0200 Subject: [PATCH 3/5] Synchronize multispecies phi refinement --- NEO-2-QL/propagator.f90 | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/NEO-2-QL/propagator.f90 b/NEO-2-QL/propagator.f90 index 3956400f..dcedc490 100644 --- a/NEO-2-QL/propagator.f90 +++ b/NEO-2-QL/propagator.f90 @@ -1589,9 +1589,13 @@ SUBROUTINE ripple_solver_int(ierr) !! Modification by Andreas F. Martitsch (27.07.2015) ! Switch for ripple_solver version USE ntv_mod, ONLY : isw_ripple_solver + USE mpi, ONLY : MPI_ALLREDUCE, MPI_COMM_WORLD, MPI_INTEGER, MPI_MAX, & + MPI_SUCCESS !! End Modification by Andreas F. Martitsch (27.07.2015) INTEGER :: ierr + INTEGER :: mpi_ierr + INTEGER, DIMENSION(2) :: local_status, global_status INTERFACE ripple_solver SUBROUTINE ripple_solver( & @@ -1672,6 +1676,26 @@ END SUBROUTINE ripple_solver_ArnoldiO2 ENDIF !! End Modification by Andreas F. Martitsch (27.07.2015) + ! In a multispecies run each MPI rank advances one species. A crossing can + ! request phi refinement on only one species, but every rank must rebuild + ! the same propagator grid before later fixed-size collectives. Promote a + ! recoverable or fatal status to all ranks, with fatal taking precedence. + local_status = 0 + IF (ierr .EQ. 3) local_status(1) = 1 + IF (ierr .NE. 0 .AND. ierr .NE. 3) local_status(2) = 1 + CALL MPI_ALLREDUCE(local_status,global_status,SIZE(local_status), & + MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,mpi_ierr) + IF (mpi_ierr .NE. MPI_SUCCESS) THEN + ierr = 1 + RETURN + ELSEIF (global_status(2) .NE. 0) THEN + ierr = 1 + ELSEIF (global_status(1) .NE. 0) THEN + ierr = 3 + ELSE + ierr = 0 + END IF + ! 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 From 24b5c72bb04d84dca05393f3aa5e59f646d6d879 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 21 Jul 2026 12:15:15 +0200 Subject: [PATCH 4/5] Synchronize refinement before solver collectives --- NEO-2-QL/propagator.f90 | 24 ------------------- NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 | 12 +++++++++- NEO-2-QL/ripple_solver_axi_test.f90 | 12 +++++++++- 3 files changed, 22 insertions(+), 26 deletions(-) diff --git a/NEO-2-QL/propagator.f90 b/NEO-2-QL/propagator.f90 index dcedc490..3956400f 100644 --- a/NEO-2-QL/propagator.f90 +++ b/NEO-2-QL/propagator.f90 @@ -1589,13 +1589,9 @@ SUBROUTINE ripple_solver_int(ierr) !! Modification by Andreas F. Martitsch (27.07.2015) ! Switch for ripple_solver version USE ntv_mod, ONLY : isw_ripple_solver - USE mpi, ONLY : MPI_ALLREDUCE, MPI_COMM_WORLD, MPI_INTEGER, MPI_MAX, & - MPI_SUCCESS !! End Modification by Andreas F. Martitsch (27.07.2015) INTEGER :: ierr - INTEGER :: mpi_ierr - INTEGER, DIMENSION(2) :: local_status, global_status INTERFACE ripple_solver SUBROUTINE ripple_solver( & @@ -1676,26 +1672,6 @@ END SUBROUTINE ripple_solver_ArnoldiO2 ENDIF !! End Modification by Andreas F. Martitsch (27.07.2015) - ! In a multispecies run each MPI rank advances one species. A crossing can - ! request phi refinement on only one species, but every rank must rebuild - ! the same propagator grid before later fixed-size collectives. Promote a - ! recoverable or fatal status to all ranks, with fatal taking precedence. - local_status = 0 - IF (ierr .EQ. 3) local_status(1) = 1 - IF (ierr .NE. 0 .AND. ierr .NE. 3) local_status(2) = 1 - CALL MPI_ALLREDUCE(local_status,global_status,SIZE(local_status), & - MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,mpi_ierr) - IF (mpi_ierr .NE. MPI_SUCCESS) THEN - ierr = 1 - RETURN - ELSEIF (global_status(2) .NE. 0) THEN - ierr = 1 - ELSEIF (global_status(1) .NE. 0) THEN - ierr = 3 - ELSE - ierr = 0 - END IF - ! 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 diff --git a/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 b/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 index d50cf804..d073acd1 100644 --- a/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 +++ b/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 diff --git a/NEO-2-QL/ripple_solver_axi_test.f90 b/NEO-2-QL/ripple_solver_axi_test.f90 index 2bdad1c3..49da592c 100644 --- a/NEO-2-QL/ripple_solver_axi_test.f90 +++ b/NEO-2-QL/ripple_solver_axi_test.f90 @@ -86,6 +86,8 @@ SUBROUTINE ripple_solver( & nz_regper,irow_regper,icol_regper,amat_regper, & deallocate_ntv_eqmat use mpiprovider_module, only : mpro + 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 +278,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 @@ -800,7 +803,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 From 02291d81d14255ebbdca6dbecb641508c6baf3fd Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 21 Jul 2026 12:56:59 +0200 Subject: [PATCH 5/5] Synchronize multispecies integral iterations --- NEO-2-QL/CMakeSources.in | 1 + NEO-2-QL/multispecies_iteration_sync_mod.f90 | 35 ++++++++++++++++ NEO-2-QL/ripple_solver_axi_test.f90 | 20 +++++---- TEST/CMakeLists.txt | 17 ++++++++ TEST/test_multispecies_iteration_sync.f90 | 43 ++++++++++++++++++++ 5 files changed, 108 insertions(+), 8 deletions(-) create mode 100644 NEO-2-QL/multispecies_iteration_sync_mod.f90 create mode 100644 TEST/test_multispecies_iteration_sync.f90 diff --git a/NEO-2-QL/CMakeSources.in b/NEO-2-QL/CMakeSources.in index 6d47d9d8..6a67daa3 100644 --- a/NEO-2-QL/CMakeSources.in +++ b/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 diff --git a/NEO-2-QL/multispecies_iteration_sync_mod.f90 b/NEO-2-QL/multispecies_iteration_sync_mod.f90 new file mode 100644 index 00000000..3add1687 --- /dev/null +++ b/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 diff --git a/NEO-2-QL/ripple_solver_axi_test.f90 b/NEO-2-QL/ripple_solver_axi_test.f90 index 49da592c..447cdf74 100644 --- a/NEO-2-QL/ripple_solver_axi_test.f90 +++ b/NEO-2-QL/ripple_solver_axi_test.f90 @@ -86,6 +86,7 @@ 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 @@ -289,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) !*************************** @@ -3136,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 @@ -3192,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 diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 0ae0e34d..494c4d35 100644 --- a/TEST/CMakeLists.txt +++ b/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} + $ + ${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) diff --git a/TEST/test_multispecies_iteration_sync.f90 b/TEST/test_multispecies_iteration_sync.f90 new file mode 100644 index 00000000..154c4e11 --- /dev/null +++ b/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