module test_parallel_map_m use pfunit use MPI implicit none contains @test(npes = [4]) subroutine parallel_map_circular(this) !! Test of parallel map for circular geometry use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP use helpers_mesh_sample_m, only : mesh_sample_circular use helpers_m, only : almost_equal, digit_sum_csr use helpers_grillix_m, only : get_mat_entry_cartbased use parallel_map_m, only : parallel_map_t use constants_m, only : two_pi use csrmat_m, only : csrmat_t use equilibrium_m, only : equilibrium_t use circular_equilibrium_m, only : circular_t use mesh_cart_m, only : mesh_cart_t class (MpiTestMethod), intent(inout) :: this integer :: comm_world type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi type(mesh_cart_t) :: mesh_cano type(mesh_cart_t) :: mesh_stag type(parallel_map_t) :: map real(GP), parameter :: trace_prec = 10E-7_GP integer :: rank, ierr real(GP) :: phi_cano, phi_stag, dphi, x, y, dpar, vol, dpar_full logical :: fine integer :: l, iret, isum, jsum real(GP) :: rret, valsum, sum_npoints comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_parallel_map_circular ------------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi_cano = dphi * comm_handler%get_plane() phi_stag = phi_cano + dphi / 2.0_GP ! Create samplemeshes ------------------------------------------------------------------------- call mesh_sample_circular(equi, phi_cano, mesh_cano) call mesh_sample_circular(equi, phi_stag, mesh_stag) ! Create map ---------------------------------------------------------------------------------- call map%set_parameters(comm_handler, & intorder_in=3, & xorder_in=1, & use_fixed_stencil_in=.false., & use_gauss_quadrature_in=.false.) call map%build(comm_handler, equi, mesh_cano, mesh_stag) ! Test meta information of map ---------------------------------------------------------------- iret = map%get_nplanes() @assertEqual(iret, 4) iret = map%get_intorder() @assertEqual(iret, 3) iret = map%get_xorder() @assertEqual(iret, 1) @assertTrue(.not.map%is_fixed_stencil()) @assertTrue(.not.map%is_gauss_quadrature()) rret = map%get_dphi() fine = almost_equal(rret, dphi, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) ! Check phi and phi_stag rret = map%get_phi_cano() fine = almost_equal(rret, phi_cano, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = map%get_phi_stag() fine = almost_equal(rret, phi_stag, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) ! Test parallel distances and flux box volumes vs. analytic result ---------------------------- do l = 1, mesh_cano%get_n_points() x = mesh_cano%get_x(l) y = mesh_cano%get_y(l) dpar = dpar_analytic(x, y, dphi / 2.0_GP) dpar_full = dpar_analytic(x, y, dphi) vol = dphi * mesh_cano%get_spacing_f()**2 rret = map%dpar_cano_stag_fwd(l) fine = almost_equal(rret, dpar, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%dpar_stag_cano_fwd(l) fine = almost_equal(rret, dpar, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%dpar_cano_stag_bwd(l) fine = almost_equal(rret, dpar, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%dpar_stag_cano_bwd(l) fine = almost_equal(rret, dpar, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%dpar_cano_cano_fwd(l) fine = almost_equal(rret, dpar_full, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%dpar_stag_stag_fwd(l) fine = almost_equal(rret, dpar_full, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%dpar_cano_cano_bwd(l) fine = almost_equal(rret, dpar_full, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%dpar_stag_stag_bwd(l) fine = almost_equal(rret, dpar_full, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%fluxbox_vol_cano(l) fine = almost_equal(rret, vol, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) rret = map%fluxbox_vol_stag(l) fine = almost_equal(rret, vol, rtol=trace_prec, atol=0.0_GP) @assertTrue(fine) enddo ! Test nnz and digit sum of map matrices ------------------------------------------------------ sum_npoints = 1.0_GP*mesh_cano%get_n_points() @assertEqual(map%map_stag_cano_fwd%nnz, 66904) call digit_sum_csr(map%map_stag_cano_fwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_stag_cano_bwd%nnz, 66904) call digit_sum_csr(map%map_stag_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_cano_stag_fwd%nnz, 66904) call digit_sum_csr(map%map_cano_stag_fwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_cano_stag_bwd%nnz, 66904) call digit_sum_csr(map%map_cano_stag_bwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%grad_stag_cano_fwd%nnz, 87804) call digit_sum_csr(map%grad_stag_cano_fwd, isum, jsum, valsum) fine = almost_equal(valsum, 3118.55513413472_GP , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%grad_stag_cano_bwd%nnz, 87804) call digit_sum_csr(map%grad_stag_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, 3118.55513413472_GP , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%pdiv_cano_stag_fwd%nnz, 87804) call digit_sum_csr(map%pdiv_cano_stag_fwd, isum, jsum, valsum) fine = almost_equal(valsum, -3118.55513413472_GP , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%pdiv_cano_stag_bwd%nnz, 87804) call digit_sum_csr(map%pdiv_cano_stag_bwd, isum, jsum, valsum) fine = almost_equal(valsum, -3118.55513413472_GP , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_cano_cano_fwd%nnz, 66928) call digit_sum_csr(map%map_cano_cano_fwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_cano_cano_fwd2%nnz, 66436) call digit_sum_csr(map%map_cano_cano_fwd2, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_cano_cano_bwd%nnz, 66928) call digit_sum_csr(map%map_cano_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_cano_cano_bwd2%nnz, 66436) call digit_sum_csr(map%map_cano_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_stag_stag_fwd%nnz, 66928) call digit_sum_csr(map%map_stag_stag_fwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) @assertEqual(map%map_stag_stag_bwd%nnz, 66928) call digit_sum_csr(map%map_stag_stag_bwd, isum, jsum, valsum) fine = almost_equal(valsum, sum_npoints , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) ! Test specific entries of map matrices (mesh order independent) ------------------------------ rret = get_mat_entry_cartbased(map%map_stag_cano_fwd, mesh_cano, & it=8, jt= 54, ig=8, jg=47) fine = almost_equal(rret, 0.809198339733931_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_stag_cano_bwd, mesh_cano, & it=8, jt=54, ig=9, jg=62) fine = almost_equal(rret, 0.118575315235283_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_cano_stag_fwd, mesh_cano, & it=8, jt=54, ig=8, jg=47) fine = almost_equal(rret, 0.809198339733931_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_cano_stag_bwd, mesh_cano, & it=8, jt=54, ig=9, jg=62) fine = almost_equal(rret, 0.118575315235283_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%grad_stag_cano_fwd, mesh_cano, & it=8, jt=54, ig=8, jg=47) fine = almost_equal(rret, 0.409722719573700_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%grad_stag_cano_bwd, mesh_cano, & it=8, jt=54, ig=9, jg=62) fine = almost_equal(rret, 9.250555723316595E-002_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%pdiv_cano_stag_fwd, mesh_cano, & it=8, jt=54, ig=8, jg=47) fine = almost_equal(rret, -0.410850057436891_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%pdiv_cano_stag_bwd, mesh_cano, & it=8, jt=54, ig=9, jg=62) fine = almost_equal(rret, -8.425052394518580E-002_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_cano_cano_fwd, mesh_cano, & it=8, jt=54, ig=10, jg=40) fine = almost_equal(rret, 0.351134771513594_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_cano_cano_fwd2, mesh_cano, & it=8, jt=54, ig=16, jg=27) fine = almost_equal(rret, 0.348793860682569_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_cano_cano_bwd, mesh_cano, & it=8, jt=54, ig=11, jg=67) fine = almost_equal(rret, 3.091610740407465E-002_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_cano_cano_bwd2, mesh_cano, & it=8, jt=54, ig=19, jg=80) fine = almost_equal(rret, 0.896898088728458_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_stag_stag_fwd, mesh_cano, & it=8, jt=54, ig=10, jg=40) fine = almost_equal(rret, 0.351134771513594_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = get_mat_entry_cartbased(map%map_stag_stag_bwd, mesh_cano, & it=8, jt=54, ig=11, jg=67) fine = almost_equal(rret, 3.091610740407465E-002_GP, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) if (rank == 0) then write(*,'(A80)')'test_parallel_map_circular complete ---------------------------------------' write(*,*)'' write(*,*)'' endif contains real(GP) function dpar_analytic(x, y, dphi) !! Returns parallel distance based on analytic calculations real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: dphi !! Toroidal distance real(GP) :: rho, qval, phi select type(equi) type is(circular_t) phi = 0.0_GP ! Circular is axisymmetric rho = equi%rho(x, y, phi) qval = equi%qval(rho) end select dpar_analytic = sqrt(1.0_GP + rho**2/qval**2) * dphi end end subroutine @test(npes = [4]) subroutine parallel_map_cerfons(this) !! test parallel_map for cerfons equilibrium use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP use helpers_mesh_sample_m, only : mesh_sample_cerfons use helpers_m, only : almost_equal, digit_sum_csr use helpers_grillix_m, only : get_mat_entry_cartbased use parallel_map_m, only : parallel_map_t use constants_m, only : two_pi use equilibrium_m, only : equilibrium_t use circular_equilibrium_m, only : circular_t use mesh_cart_m, only : mesh_cart_t class (MpiTestMethod), intent(inout) :: this integer :: comm_world type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi type(mesh_cart_t) :: mesh_cano type(mesh_cart_t) :: mesh_stag type(parallel_map_t) :: map logical :: fine integer :: l, iret, isum, jsum, rank, ierr real(GP) :: rret, phi_cano, phi_stag, dphi, & sum_dparfwd, sum_dparbwd, sum_vol, valsum, sum_dparfwd_full, sum_dparbwd_full, & sum_dparfwd_stag, sum_dparbwd_stag, sum_vol_stag, sum_dparfwd_full_stag, sum_dparbwd_full_stag comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_parallel_map_cerfons -------------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi_cano = dphi * comm_handler%get_plane() phi_stag = phi_cano + dphi / 2.0_GP ! Create samplemesh --------------------------------------------------------------------------- call mesh_sample_cerfons(equi, phi_cano, mesh_cano) call mesh_sample_cerfons(equi, phi_stag, mesh_stag) ! Create map ---------------------------------------------------------------------------------- call map%set_parameters(comm_handler, & intorder_in = 3, & xorder_in = 1, & use_fixed_stencil_in = .false., & use_gauss_quadrature_in = .false. ) call map%build(comm_handler, equi, mesh_cano, mesh_stag) ! Test meta information of map ---------------------------------------------------------------- iret = map%get_nplanes() @assertEqual(iret, 4) iret = map%get_intorder() @assertEqual(iret, 3) iret = map%get_xorder() @assertEqual(iret, 1) @assertTrue(.not.map%is_fixed_stencil()) @assertTrue(.not.map%is_gauss_quadrature()) ! Check metadata rret = map%get_dphi() fine = almost_equal(rret, dphi, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = map%get_phi_cano() fine = almost_equal(rret, phi_cano, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rret = map%get_phi_stag() fine = almost_equal(rret, phi_stag, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) ! Test Digit sum if dpar and vol -------------------------------------------------------------- sum_dparfwd = 0.0_GP sum_dparbwd = 0.0_GP sum_dparfwd_full = 0.0_GP sum_dparbwd_full = 0.0_GP sum_vol = 0.0_GP sum_dparfwd_stag = 0.0_GP sum_dparbwd_stag = 0.0_GP sum_dparfwd_full_stag = 0.0_GP sum_dparbwd_full_stag = 0.0_GP sum_vol_stag = 0.0_GP do l = 1, mesh_cano%get_n_points() sum_dparfwd = sum_dparfwd + map%dpar_cano_stag_fwd(l) sum_dparbwd = sum_dparbwd + map%dpar_cano_stag_bwd(l) sum_dparfwd_full = sum_dparfwd_full + map%dpar_cano_cano_fwd(l) sum_dparbwd_full = sum_dparbwd_full + map%dpar_cano_cano_bwd(l) sum_vol = sum_vol + map%fluxbox_vol_cano(l) sum_dparfwd_stag = sum_dparfwd_stag + map%dpar_stag_cano_fwd(l) sum_dparbwd_stag = sum_dparbwd_stag + map%dpar_stag_cano_bwd(l) sum_dparfwd_full_stag = sum_dparfwd_full_stag + map%dpar_stag_stag_fwd(l) sum_dparbwd_full_stag = sum_dparbwd_full_stag + map%dpar_stag_stag_bwd(l) sum_vol_stag = sum_vol_stag + map%fluxbox_vol_stag(l) enddo fine = almost_equal(sum_dparfwd, 6694.81611206796_GP , rtol=0.0_GP, atol=1.0E-4_GP) @assertTrue(fine) fine = almost_equal(sum_dparbwd, 6816.04201374531_GP , rtol=0.0_GP, atol=1.0E-4_GP) @assertTrue(fine) fine = almost_equal(sum_dparfwd_full, 13308.0174895849_GP , rtol=0.0_GP, atol=1.0E-4_GP) @assertTrue(fine) fine = almost_equal(sum_dparbwd_full, 13692.1315956905_GP , rtol=0.0_GP, atol=1.0E-4_GP) @assertTrue(fine) fine = almost_equal(sum_vol, 1.31647650966758_GP , rtol=0.0_GP, atol=1.0E-9_GP) @assertTrue(fine) ! Staggered metadata is the same due to axisymmetry fine = almost_equal(sum_dparfwd_stag, sum_dparfwd , rtol=0.0_GP, atol=1.0E-12_GP) @assertTrue(fine) fine = almost_equal(sum_dparbwd_stag, sum_dparbwd , rtol=0.0_GP, atol=1.0E-12_GP) @assertTrue(fine) fine = almost_equal(sum_dparfwd_full_stag, sum_dparfwd_full , rtol=0.0_GP, atol=1.0E-12_GP) @assertTrue(fine) fine = almost_equal(sum_dparbwd_full_stag, sum_dparbwd_full , rtol=0.0_GP, atol=1.0E-12_GP) @assertTrue(fine) fine = almost_equal(sum_vol_stag, sum_vol , rtol=0.0_GP, atol=1.0E-12_GP) @assertTrue(fine) ! Test nnz and digit sum of map matrices ------------------------------------------------------ @assertEqual(map%map_stag_cano_fwd%nnz, 145412) call digit_sum_csr(map%map_stag_cano_fwd, isum, jsum, valsum) fine = almost_equal(valsum, 9995.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%map_stag_cano_bwd%nnz, 145032) call digit_sum_csr(map%map_stag_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, 9927.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%map_cano_stag_fwd%nnz, 145412) call digit_sum_csr(map%map_stag_cano_fwd, isum, jsum, valsum) fine = almost_equal(valsum, 9995.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%map_cano_stag_bwd%nnz, 145032) call digit_sum_csr(map%map_stag_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, 9927.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%grad_stag_cano_fwd%nnz, 196917) call digit_sum_csr(map%grad_stag_cano_fwd, isum, jsum, valsum) fine = almost_equal(valsum, 8101.82241142591_GP, rtol=0.0_GP, atol=1.0E-6_GP) @assertTrue(fine) @assertEqual(map%grad_stag_cano_bwd%nnz, 196479) call digit_sum_csr(map%grad_stag_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, 8117.60998935224_GP, rtol=0.0_GP, atol=1.0E-6_GP) @assertTrue(fine) @assertEqual(map%pdiv_cano_stag_fwd%nnz, 196479) call digit_sum_csr(map%pdiv_cano_stag_fwd, isum, jsum, valsum) fine = almost_equal(valsum, -8081.70498408458_GP, rtol=0.0_GP, atol=1.0E-6_GP) @assertTrue(fine) @assertEqual(map%pdiv_cano_stag_bwd%nnz, 196917) call digit_sum_csr(map%pdiv_cano_stag_bwd, isum, jsum, valsum) fine = almost_equal(valsum, -8176.19107533937_GP, rtol=0.0_GP, atol=1.0E-6_GP) @assertTrue(fine) @assertEqual(map%map_cano_cano_fwd%nnz, 143747) call digit_sum_csr(map%map_cano_cano_fwd, isum, jsum, valsum) fine = almost_equal(valsum, 9860.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%map_cano_cano_fwd2%nnz, 139604) call digit_sum_csr(map%map_cano_cano_fwd2, isum, jsum, valsum) fine = almost_equal(valsum, 9494.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%map_cano_cano_bwd%nnz, 142906) call digit_sum_csr(map%map_cano_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, 9751.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%map_cano_cano_bwd2%nnz, 138160) call digit_sum_csr(map%map_cano_cano_bwd2, isum, jsum, valsum) fine = almost_equal(valsum, 9319.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%map_stag_stag_fwd%nnz, 143747) call digit_sum_csr(map%map_cano_cano_fwd, isum, jsum, valsum) fine = almost_equal(valsum, 9860.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) @assertEqual(map%map_stag_stag_bwd%nnz, 142906) call digit_sum_csr(map%map_cano_cano_bwd, isum, jsum, valsum) fine = almost_equal(valsum, 9751.0_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) ! Note: Matrices itself are also well tested within PARALLAX if (rank == 0) then write(*,'(A80)')'test_parallel_map_cerfons complete ---------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine end module