module test_parallel_target_flux_m !! Tests of parallel target flux module use pfunit use netcdf use MPI use constants_m, only : TWO_PI use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP use equilibrium_m, only : equilibrium_t use equilibrium_factory_m, only : create_equilibrium use equilibrium_storage_m, only : equilibrium_storage_t use circular_equilibrium_m, only : circular_t use cerfons_equilibrium_m, only : cerfons_t use mesh_cart_m, only : mesh_cart_t use multigrid_m, only : multigrid_t use parallel_map_m, only : parallel_map_t use penalisation_m, only : penalisation_t use helpers_mesh_sample_m, only : mesh_sample_circular, mesh_sample_cerfons use helpers_m, only : almost_equal use screen_io_m, only : get_stdout use error_handling_grillix_m, only : handle_error_netcdf use parallel_target_flux_m, only : parallel_target_flux_t implicit none integer, private, parameter :: HERMITE_PF = 465 integer, private, parameter :: VIA_TRACE = 811 contains @test(npes = [4]) subroutine parallel_target_flux_circular(this) !! Tests parallel target flux in CIRCULAR geometry class (MpiTestMethod), intent(inout) :: this integer :: comm_world, rank, ierr type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi type(mesh_cart_t) :: mesh, mesh_stag type(equilibrium_storage_t) :: equi_st type(multigrid_t) :: multigrid type(parallel_map_t) :: map type(penalisation_t) :: penalisation type(parallel_target_flux_t) :: partargetflux logical :: fine real(GP) :: phi, dphi, phi_stag, x, y, rhon, flux_fwd, flux_bwd integer :: l real(GP), dimension(:), allocatable :: g comm_world = this%getMpiCommunicator() call MPI_COMM_RANK(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_parallel_target_flux_circular ------------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi = dphi * comm_handler%get_plane() phi_stag = phi + dphi / 2.0_GP ! Create samplemesh --------------------------------------------------------------------------- call mesh_sample_circular(equi, phi, mesh, multigrid) call mesh_sample_circular(equi, phi_stag, mesh_stag) call equi_st%fill_storage(equi, mesh) ! 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, mesh_stag) ! Create penalisation with hermite ------------------------------------------------------------ call penalisation%set_parameters(dphi, & pen_method_in = VIA_TRACE, & penfuns_type_in = HERMITE_PF, & hermite_order_in = 3, & charfun_parwidth_in = 1.2_GP, & charfun_radlimwidth_in = 5.0E-2_GP, & dirindfun_parwidth_in = 2.0_GP, & dphi_max_in = 4.0_GP*dphi, & charfun_at_target_in = 0.5_GP, & max_step_size_in = 4.0_GP*dphi, & rho_min_in = 0.0_GP) call penalisation%build(comm_handler, equi, mesh, multigrid) allocate(g(mesh%get_n_points())) ! Set up example fields do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rhon = ( equi%rho(x,y,phi) - equi%rhomin ) / ( equi%rhomax - equi%rhomin ) g(l) = ( sin(0.6_GP*TWO_PI*rhon + 0.1_GP)*cos(2.1*phi - 0.3_GP) + 1.2_GP ) & * ( 3.1_GP*sin(1.9_GP*TWO_PI*rhon - 2.5_GP)*sin(0.3_GP*phi+1.2_GP) ) enddo ! Create target markers on canonical grid call partargetflux%build(comm_handler, mesh, map, penalisation, staggered=.false., & pen_threshold=0.5_GP) ! Test number of marked indices @assertEqual(partargetflux%n_inds_fwd, 84) @assertEqual(partargetflux%n_inds_bwd, 84) ! Test target flux computation on canonical grid call partargetflux%compute(comm_handler, mesh, penalisation, equi_st, g, & flux_fwd, flux_bwd, mode='DEFAULT') fine = almost_equal(flux_fwd, 2.214874854978443E-002_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) fine = almost_equal(flux_bwd, -2.214874854978456E-002_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) call partargetflux%compute(comm_handler, mesh, penalisation, equi_st, g, & flux_fwd, flux_bwd, mode='WEIGHTED') fine = almost_equal(flux_fwd, 2.268702472864539E-002_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) fine = almost_equal(flux_bwd, -2.268702472864539E-002_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) if (rank == 0) then write(*,'(A80)')'test_parallel_target_flux_circular complete ---------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine @test(npes = [4]) subroutine parallel_target_flux_cerfons(this) !! Tests parallel target flux in CERFONS geometry class (MpiTestMethod), intent(inout) :: this integer :: comm_world, rank, ierr type(comm_handler_t) ::comm_handler class(equilibrium_t), allocatable :: equi type(mesh_cart_t) :: mesh, mesh_stag type(equilibrium_storage_t) :: equi_st type(multigrid_t) :: multigrid type(parallel_map_t) :: map type(penalisation_t) :: penalisation type(parallel_target_flux_t) :: partargetflux logical :: fine real(GP) :: phi, dphi, phi_stag, x, y, rhon, flux_fwd, flux_bwd integer :: l real(GP), dimension(:), allocatable :: g comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_parallel_target_flux_cerfons -------------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi = dphi * comm_handler%get_plane() phi_stag = phi + dphi / 2.0_GP ! Create samplemesh --------------------------------------------------------------------------- call mesh_sample_cerfons(equi, phi, mesh, multigrid) call mesh_sample_cerfons(equi, phi_stag, mesh_stag) call equi_st%fill_storage(equi, mesh) ! 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, mesh_stag) ! Create penalisation with hermite ------------------------------------------------------------ call penalisation%set_parameters(dphi, & pen_method_in = VIA_TRACE, & penfuns_type_in = HERMITE_PF, & hermite_order_in = 3, & charfun_parwidth_in = 1.0_GP, & charfun_radlimwidth_in = 0.0_GP, & dirindfun_parwidth_in = 0.8_GP, & dphi_max_in = 4.0_GP*dphi, & charfun_at_target_in = 0.5_GP, & max_step_size_in = 4.0_GP*dphi, & rho_min_in = 0.0_GP) call penalisation%build(comm_handler, equi, mesh, multigrid) allocate(g(mesh%get_n_points())) ! Set up example fields do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rhon = ( equi%rho(x,y,phi) - equi%rhomin ) / ( equi%rhomax - equi%rhomin ) g(l) = ( sin(0.6_GP*TWO_PI*rhon + 0.1_GP)*cos(2.1*phi - 0.3_GP) + 1.2_GP ) & * ( 3.1_GP*sin(1.9_GP*TWO_PI*rhon - 2.5_GP)*sin(0.3_GP*phi+1.2_GP) ) enddo ! Create target markers on canonical grid call partargetflux%build(comm_handler, mesh, map, penalisation, staggered=.false., & pen_threshold=0.5_GP) ! Test number of marked indices @assertEqual(partargetflux%n_inds_fwd, 138) @assertEqual(partargetflux%n_inds_bwd, 122) ! Test target flux computation on canonical grid call partargetflux%compute(comm_handler, mesh, penalisation, equi_st, g, & flux_fwd, flux_bwd, mode = 'DEFAULT') fine = almost_equal(flux_fwd, 9.815055630709653E-002_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) fine = almost_equal(flux_bwd, -7.826238455294574E-002_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) call partargetflux%compute(comm_handler, mesh, penalisation, equi_st, g, & flux_fwd, flux_bwd, mode = 'WEIGHTED') fine = almost_equal(flux_fwd, 9.768507970778310E-002_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) fine = almost_equal(flux_bwd, -7.774840359607528E-002_GP, rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) if (rank == 0) then write(*,'(A80)')'test_parallel_target_flux_cerfons complete ----------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine end module