module test_penalisation_m !! Tests for most of the penalisation module use pfunit use grillix_build_info_m, only : git_repo_path 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, DOMMASCHK 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 penalisation_m, only : penalisation_t use helpers_mesh_sample_m, only : mesh_sample_circular, mesh_sample_cerfons use helpers_m, only : almost_equal implicit none integer, private, parameter :: TANH_PF = 464 integer, private, parameter :: HERMITE_PF = 465 integer, private, parameter :: VIA_TRACE = 811 integer, private, parameter :: VIA_IN_VESSEL = 812 integer, private, parameter :: VIA_NONE = 813 integer, private, parameter :: VIA_STABLE_TRACE = 815 contains @test(npes = [4]) subroutine penalisation_circular(this) !! Tests penalisation 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 type(multigrid_t) :: multigrid type(penalisation_t) :: pentanh, pentanh_shifted type(penalisation_t) :: penhermite, penhermite_shifted logical :: fine integer :: ki, iret, it, jt, lt real(GP) :: phi, dphi, rret comm_world = this%getMpiCommunicator() call MPI_COMM_RANK(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_penalisation_circular ------------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi = dphi * comm_handler%get_plane() ! Create samplemesh --------------------------------------------------------------------------- call mesh_sample_circular(equi, phi, mesh, multigrid) ! Create penalisation with TANH---------------------------------------------------------------- call pentanh%set_parameters(dphi, & pen_method_in = VIA_TRACE, & penfuns_type_in = TANH_PF, & hermite_order_in = 3, & charfun_parwidth_in = 0.8_GP, & charfun_radlimwidth_in = 2.0E-2_GP, & dirindfun_parwidth_in = 0.6_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 pentanh%build(comm_handler, equi, mesh, multigrid) ! Create additional penalisation with shifted charfun values call pentanh_shifted%set_parameters(dphi, & pen_method_in = VIA_TRACE, & penfuns_type_in = TANH_PF, & hermite_order_in = 3, & charfun_parwidth_in = 0.8_GP, & charfun_radlimwidth_in = 2.0E-2_GP, & dirindfun_parwidth_in = 0.6_GP , & dphi_max_in = 4.0_GP*dphi, & charfun_at_target_in = 0.3_GP, & max_step_size_in = 4.0_GP*dphi, & rho_min_in = 0.0_GP) call pentanh_shifted%build(comm_handler, equi, mesh, multigrid) ! Create penalisation with hermite ------------------------------------------------------------ call penhermite%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 penhermite%build(comm_handler, equi, mesh, multigrid) ! Create additional penalisation with shifted charfun values call penhermite_shifted%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.1_GP, & max_step_size_in = 4.0_GP*dphi, & rho_min_in = 0.0_GP) call penhermite_shifted%build(comm_handler, equi, mesh, multigrid) ! Test meta information of pentanh ------------------------------------------------------------ iret = pentanh%get_penfuns_type() @assertEqual(iret, TANH_PF) iret = pentanh%get_hermite_order() @assertEqual(iret, 3) rret = pentanh%get_charfun_parwidth() fine = almost_equal(rret, 0.8_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_charfun_radlimwidth() fine = almost_equal(rret, 2.0E-2_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_dirindfun_parwidth() fine = almost_equal(rret, 0.6_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_dphi() fine = almost_equal(rret, dphi, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_dphi_max() fine = almost_equal(rret, 8.0_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_charfun_at_target() fine = almost_equal(rret, 0.5_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh_shifted%get_charfun_at_target() fine = almost_equal(rret, 0.3_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! Test meta information of penhermite --------------------------------------------------------- iret = penhermite%get_penfuns_type() @assertEqual(iret, HERMITE_PF) iret = penhermite%get_hermite_order() @assertEqual(iret, 3) rret = penhermite%get_charfun_parwidth() fine = almost_equal(rret, 1.2_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_charfun_radlimwidth() fine = almost_equal(rret, 5.0E-2_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_dirindfun_parwidth() fine = almost_equal(rret, 2.0_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_dphi() fine = almost_equal(rret, dphi, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_dphi_max() fine = almost_equal(rret, 4.0_GP*dphi, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_charfun_at_target() fine = almost_equal(rret, 0.5_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite_shifted%get_charfun_at_target() fine = almost_equal(rret, 0.1_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! Test values for a specific point for pentanh and penhermite --------------------------------- it = 8 jt = 54 lt = mesh%cart_to_mesh_index(it,jt) fine = .false. do ki = 1, mesh%get_n_points_inner() if (mesh%inner_indices(ki) == lt) then rret = pentanh%get_charfun_val(ki) fine = almost_equal(rret, 0.979716115268211_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) rret = pentanh_shifted%get_charfun_val(ki) fine = almost_equal(rret, 0.953917246613992_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) rret = pentanh%get_dirindfun_val(ki) fine = almost_equal(rret, 0.796271722032819_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) rret = penhermite%get_charfun_val(ki) fine = almost_equal(rret, 0.343106979541614_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) rret = penhermite_shifted%get_charfun_val(ki) fine = almost_equal(rret, 0.343106979541614_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) rret = penhermite%get_dirindfun_val(ki) fine = almost_equal(rret, 0.959354522571291_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) endif enddo ! Check if point was found at all @assertTrue(fine) ! Test pen index storage sizes @assertEqual(pentanh%n_p_inds, 1458) @assertEqual(pentanh%n_pb_inds, 324) @assertEqual(penhermite%n_p_inds, 239) @assertEqual(penhermite%n_pb_inds, 65) ! Compute checksums for penalisation indices iret = pen_ind_checksum(pentanh%p_inds) @assertEqual(iret, 3579724) iret = pen_ind_checksum(pentanh%pb_inds) @assertEqual(iret, 797984) iret = pen_ind_checksum(penhermite%p_inds) @assertEqual(iret, 390925) iret = pen_ind_checksum(penhermite%pb_inds) @assertEqual(iret, 108284) ! Check some individual indices iret = pentanh%p_inds(29) @assertEqual(iret, 68) iret = pentanh%pb_inds(7) @assertEqual(iret, 220) iret = penhermite%p_inds(29) @assertEqual(iret, 1102) iret = penhermite%pb_inds(7) @assertEqual(iret, 318) if (rank == 0) then write(*,'(A80)')'test_penalisation_circular complete ---------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine @test(npes = [4]) subroutine penalisation_cerfons(this) !! Tests penalisation 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 type(multigrid_t) :: multigrid type(penalisation_t) :: pentanh, pentanh_shifted type(penalisation_t) :: penhermite, penhermite_shifted logical :: fine integer :: ki, iret, it, jt, lt real(GP) :: phi, dphi, rret comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_penalisation_cerfons -------------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi = dphi * comm_handler%get_plane() ! Create samplemesh --------------------------------------------------------------------------- call mesh_sample_cerfons(equi, phi, mesh, multigrid) ! Create penalisation with TANH---------------------------------------------------------------- call pentanh%set_parameters(dphi, & pen_method_in = VIA_TRACE, & penfuns_type_in = TANH_PF, & hermite_order_in = 3, & charfun_parwidth_in = 0.4_GP, & charfun_radlimwidth_in = 0.0_GP, & dirindfun_parwidth_in = 0.4_GP, & dphi_max_in = 2.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 pentanh%build(comm_handler, equi, mesh, multigrid) ! Create additional penalisation with shifted charfun values call pentanh_shifted%set_parameters(dphi, & pen_method_in = VIA_TRACE, & penfuns_type_in = TANH_PF, & hermite_order_in = 3, & charfun_parwidth_in = 0.4_GP, & charfun_radlimwidth_in = 0.0_GP, & dirindfun_parwidth_in = 0.4_GP, & dphi_max_in = 2.0_GP*dphi, & charfun_at_target_in = 0.3_GP, & max_step_size_in = 4.0_GP*dphi, & rho_min_in = 0.0_GP) call pentanh_shifted%build(comm_handler, equi, mesh, multigrid) ! Create penalisation with hermite ------------------------------------------------------------ call penhermite%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 = 2.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 penhermite%build(comm_handler, equi, mesh, multigrid) ! Create additional penalisation with shifted charfun values call penhermite_shifted%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 = 2.0_GP*dphi, & charfun_at_target_in = 0.3_GP, & max_step_size_in = 4.0_GP*dphi, & rho_min_in = 0.0_GP) call penhermite_shifted%build(comm_handler, equi, mesh, multigrid) ! Test meta information of pentanh ------------------------------------------------------------ iret = pentanh%get_penfuns_type() @assertEqual(iret, TANH_PF) iret = pentanh%get_hermite_order() @assertEqual(iret, 3) rret = pentanh%get_charfun_parwidth() fine = almost_equal(rret, 0.4_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_charfun_radlimwidth() fine = almost_equal(rret, 0.0_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_dirindfun_parwidth() fine = almost_equal(rret, 0.4_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_dphi() fine = almost_equal(rret, dphi, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_dphi_max() fine = almost_equal(rret, 4.0_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh%get_charfun_at_target() fine = almost_equal(rret, 0.5_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = pentanh_shifted%get_charfun_at_target() fine = almost_equal(rret, 0.3_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! Test meta information of penhermite --------------------------------------------------------- iret = penhermite%get_penfuns_type() @assertEqual(iret, HERMITE_PF) iret = penhermite%get_hermite_order() @assertEqual(iret, 3) rret = penhermite%get_charfun_parwidth() fine = almost_equal(rret, 1.0_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_charfun_radlimwidth() fine = almost_equal(rret, 0.0_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_dirindfun_parwidth() fine = almost_equal(rret, 0.8_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_dphi() fine = almost_equal(rret, dphi, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_dphi_max() fine = almost_equal(rret, 2.0_GP*dphi, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_charfun_at_target() fine = almost_equal(rret, 0.5_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite_shifted%get_charfun_at_target() fine = almost_equal(rret, 0.3_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! Test values of chi for a specific point for pentanh and penhermite -------------------------- it = 22 jt = 10 lt = mesh%cart_to_mesh_index(it,jt) fine = .false. do ki = 1, mesh%get_n_points_inner() if (mesh%inner_indices(ki) == lt) then rret = pentanh%get_charfun_val(ki) fine = almost_equal(rret, 0.4637990637803_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) rret = pentanh_shifted%get_charfun_val(ki) fine = almost_equal(rret, 0.270447083499598_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) rret = penhermite%get_charfun_val(ki) fine = almost_equal(rret, 0.4367505382371_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) rret = penhermite_shifted%get_charfun_val(ki) fine = almost_equal(rret, 0.245240794918281_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) endif enddo ! Check if point was found at all @assertTrue(fine) ! Test values of xi for a specific point for pentanh and penhermite --------------------------- it = 33 jt = 15 lt = mesh%cart_to_mesh_index(it,jt) fine = .false. do ki = 1, mesh%get_n_points_inner() if (mesh%inner_indices(ki) == lt) then rret = pentanh%get_dirindfun_val(ki) fine = almost_equal(rret, -3.139881146725E-001_GP, rtol = 0.0_GP, atol = 5.0E-6_GP) @assertTrue(fine) rret = penhermite%get_dirindfun_val(ki) fine = almost_equal(rret, -6.404270810507E-001_GP, rtol = 0.0_GP, atol = 5.0E-6_GP) @assertTrue(fine) endif enddo ! Check if point was found at all @assertTrue(fine) ! Test pen index storage sizes @assertEqual(pentanh%n_p_inds, 8604) @assertEqual(pentanh%n_pb_inds, 0) @assertEqual(penhermite%n_p_inds, 206) @assertEqual(penhermite%n_pb_inds, 31) ! Compute checksums for penalisation indices iret = pen_ind_checksum(pentanh%p_inds) @assertEqual(iret, 43743408) iret = pen_ind_checksum(pentanh%pb_inds) @assertEqual(iret, 0) iret = pen_ind_checksum(penhermite%p_inds) @assertEqual(iret, 69591) iret = pen_ind_checksum(penhermite%pb_inds) @assertEqual(iret, 13748) ! Check some individual indices iret = pentanh%p_inds(57) @assertEqual(iret, 124) iret = penhermite%p_inds(29) @assertEqual(iret, 78) iret = penhermite%pb_inds(7) @assertEqual(iret, 212) if (rank == 0) then write(*,'(A80)')'test_penalisation_cerfons complete ----------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine @test(npes = [4]) subroutine penalisation_dommaschk(this) !! Tests penalisation in DOMMASCHK geometry !! only VIA_IN_VESSEL and NONE methods class (MpiTestMethod), intent(inout) :: this integer :: comm_world, rank, ierr type(comm_handler_t) ::comm_handler class(equilibrium_t), allocatable :: equi character(len=:), allocatable :: equi_params_file type(mesh_cart_t) :: mesh type(multigrid_t) :: multigrid type(penalisation_t) :: pen_in_vessel type(penalisation_t) :: pen_none logical :: fine integer :: ki, l real(GP) :: spacing_f, phi, dphi, x, y, rref, rret comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_penalisation_dommasck ------------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi = dphi * comm_handler%get_plane() ! Create equilibrium equi_params_file = git_repo_path // & '/cicd/res/params_equi_dommaschk.in' call create_equilibrium(equi, DOMMASCHK, equi_params_file) ! Create mesh spacing_f = 1E-2 call multigrid%init(equi, phi, 2, spacing_f, 2, 2, mesh, [8,8], .false., 0) ! Create penalisation with IN_VESSEL call pen_in_vessel%set_parameters(dphi, & pen_method_in = VIA_IN_VESSEL, & penfuns_type_in = TANH_PF, & ! from here dummy parameters hermite_order_in = 3, & charfun_parwidth_in = 0.0_GP, & charfun_radlimwidth_in = 0.0_GP, & dirindfun_parwidth_in = 0.0_GP, & dphi_max_in = 0.0_GP, & charfun_at_target_in = 0.5_GP, & max_step_size_in = 4.0_GP*dphi, & rho_min_in = 0.0_GP) call pen_in_vessel%build(comm_handler, equi, mesh, multigrid) do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) x = mesh%get_x(l) y = mesh%get_y(l) if (equi%in_vessel(x, y, phi)) then rref = 0.0_GP else rref = 1.0_GP endif rret = pen_in_vessel%get_charfun_val(ki) fine = almost_equal(rref, rret, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rref = 0.0_GP rret = pen_in_vessel%get_dirindfun_val(ki) fine = almost_equal(rref, rret, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) enddo ! Test pen index storage sizes @assertEqual(pen_in_vessel%n_p_inds, 0) @assertEqual(pen_in_vessel%n_pb_inds, 0) ! Create penalisation with NONE call pen_none%set_parameters(dphi, & pen_method_in = VIA_NONE, & penfuns_type_in = TANH_PF, & ! from here dummy parameters hermite_order_in = 3, & charfun_parwidth_in = 0.0_GP, & charfun_radlimwidth_in = 0.0_GP, & dirindfun_parwidth_in = 0.0_GP, & dphi_max_in = 0.0_GP, & charfun_at_target_in = 0.5_GP, & max_step_size_in = 4.0_GP*dphi, & rho_min_in = 0.0_GP) call pen_none%build(comm_handler, equi, mesh, multigrid) do ki = 1, mesh%get_n_points_inner() rref = 0.0_GP rret = pen_none%get_charfun_val(ki) fine = almost_equal(rref, rret, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) rref = 0.0_GP rret = pen_none%get_dirindfun_val(ki) fine = almost_equal(rref, rret, rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) enddo ! Test pen index storage sizes @assertEqual(pen_none%n_p_inds, 0) @assertEqual(pen_none%n_pb_inds, 0) if (rank == 0) then write(*,'(A80)')'test_penalisation_dommaschk complete ---------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine @test(npes = [4]) subroutine penalisation_circular_stable_trace(this) !! Tests penalisation in CIRCULAR geometry via the stable trace method 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 type(multigrid_t) :: multigrid type(penalisation_t) :: penhermite logical :: fine integer :: ki, iret, it, jt, lt real(GP) :: phi, dphi, rret comm_world = this%getMpiCommunicator() call MPI_COMM_RANK(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_stable_penalisation_circular -----------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi = dphi * comm_handler%get_plane() ! Create samplemesh --------------------------------------------------------------------------- call mesh_sample_circular(equi, phi, mesh, multigrid) ! Create penalisation with hermite ------------------------------------------------------------ call penhermite%set_parameters(dphi, & pen_method_in = VIA_STABLE_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 = 0.01_GP, & rho_min_in = 0.32_GP) call penhermite%build(comm_handler, equi, mesh, multigrid) ! Test values for a specific point for pentanh and penhermite --------------------------------- it = 8 jt = 54 lt = mesh%cart_to_mesh_index(it,jt) fine = .false. do ki = 1, mesh%get_n_points_inner() if (mesh%inner_indices(ki) == lt) then rret = penhermite%get_charfun_val(ki) fine = almost_equal(rret, 1.0_GP, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) rret = penhermite%get_dirindfun_val(ki) fine = almost_equal(rret, 0.958101494873045_GP, rtol = 0.0_GP, atol = 1.0E-6_GP) @assertTrue(fine) endif enddo ! Check if point was found at all @assertTrue(fine) ! Test pen index storage sizes @assertEqual(penhermite%n_p_inds, 239) @assertEqual(penhermite%n_pb_inds, 65) ! Compute checksums for penalisation indices iret = pen_ind_checksum(penhermite%p_inds) @assertEqual(iret, 390925) iret = pen_ind_checksum(penhermite%pb_inds) @assertEqual(iret, 108284) ! Check some individual indices iret = penhermite%p_inds(29) @assertEqual(iret, 1102) iret = penhermite%pb_inds(7) @assertEqual(iret, 318) if (rank == 0) then write(*,'(A80)')'test_stable_penalisation_circular complete --------------------------------' write(*,*)'' write(*,*)'' endif end subroutine function pen_ind_checksum(u) result(cs) !! Compute sum of indices integer, dimension(:), intent(in) :: u ! Input index array integer :: cs ! Output checksum integer :: i cs = 0 !$omp parallel do private(i) reduction(+:cs) do i = 1, size(u) cs = cs + u(i) enddo !$omp end parallel do end function end module