module test_diagnostics_helpers_m use pfunit use MPI implicit none integer, private, parameter :: HERMITE_PF = 465 integer, private, parameter :: VIA_TRACE = 811 contains @test(npes = [4]) subroutine diagnostics_helpers(this) !! Test of diagnostics helper functions use precision_grillix_m, only : GP use comm_handler_m, only : comm_handler_t use helpers_mesh_sample_m, only : mesh_sample_circular use helpers_m, only : almost_equal use parallel_map_m, only : parallel_map_t use penalisation_m, only : penalisation_t ! from PARALLAX 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 use multigrid_m, only : multigrid_t ! from module to test use diagnostics_helpers_m, only : inner_product, total_volume, & interpolate_lineout implicit none 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(multigrid_t) :: multigrid type(parallel_map_t) :: map type(penalisation_t) :: penalisation logical :: fine integer :: l, n_points, int_order, ki real(GP) :: phi, dphi, phi_stag, x, y, rho, rhon, theta, res, ref, vol, vol_full real(GP), allocatable, dimension(:) :: u, v, u_inner, v_inner real(GP), dimension(5) :: x_coords, y_coords, res_lineout real(GP), dimension(5,4) :: reference_lineout comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)') & 'test_diagnostics_helpers --------------------------------------------------' 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) ! 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) ! Create some values allocate(u(mesh%get_n_points()), v(mesh%get_n_points())) allocate(u_inner(mesh%get_n_points_inner()), & v_inner(mesh%get_n_points_inner())) select type(equi) type is(circular_t) do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) theta = equi%theta(x,y) u(l) = sin(two_pi*rhon)*sin(1.2_GP*theta-0.9_GP) & * (1.0_GP + sin(phi)) v(l) = (1.5_GP*sin(2.9_GP*two_pi*rhon)*sin(5*theta-0.2_GP) + 0.5_GP) & * (1.0_GP + cos(phi)) enddo ! Fill inner_point arrays do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) u_inner(ki) = u(l) v_inner(ki) = v(l) enddo class default @assertTrue(.false.) end select ! Test total volume function res = total_volume(comm_handler, mesh, map, penalisation, & is_stag=.false., use_full_domain=.false.) vol = 1.39014218284286_GP fine = almost_equal(res, vol, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = total_volume(comm_handler, mesh, map, penalisation, & is_stag=.false., use_full_domain=.true.) vol_full = 1.97523239864754_GP fine = almost_equal(res, vol_full, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! Test inner product function ! Two input vectors ref = 1.826476260253523E-002_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, v_in=v_inner, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.false., & use_vol_avg=.true., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, v_in=v_inner, & use_abs=.false., & use_vol_avg=.true., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = 0.405352207796452_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.true., & use_vol_avg=.false., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, v_in=v_inner, & use_abs=.true., & use_vol_avg=.false., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.true., & use_vol_avg=.true., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, v_in=v_inner, & use_abs=.true., & use_vol_avg=.true., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! One input vector ref = 2.254756555463476E-002_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, & use_abs=.false., & use_vol_avg=.true., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, & use_abs=.false., & use_vol_avg=.true., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = 0.546000769866837_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, & use_abs=.true., & use_vol_avg=.false., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, & use_abs=.true., & use_vol_avg=.false., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, & use_abs=.true., & use_vol_avg=.true., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, & use_abs=.true., & use_vol_avg=.true., & use_full_domain=.false.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! With full domain enabled ! ..two input vectors ! ....full mesh input ref = -2.510160711445213E-002_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol_full res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.false., & use_vol_avg=.true., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = 0.516864889609493_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.true., & use_vol_avg=.false., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol_full res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.true., & use_vol_avg=.true., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! ....inner mesh input ref = -2.510160711445213E-002_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol_full res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.false., & use_vol_avg=.true., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = 0.516864889609493_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.true., & use_vol_avg=.false., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol_full res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, v_in=v, & use_abs=.true., & use_vol_avg=.true., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! ..one input vector ! ....full mesh input ref = 2.302368732702139E-003_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol_full res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, & use_abs=.false., & use_vol_avg=.true., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = 0.683598444762267_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, & use_abs=.true., & use_vol_avg=.false., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol_full res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u, & use_abs=.true., & use_vol_avg=.true., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! ....inner mesh input ref = 4.240319324091564E-003_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol_full res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, & use_abs=.false., & use_vol_avg=.true., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = 0.571713272540560_GP res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, & use_abs=.true., & use_vol_avg=.false., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ref = ref / vol_full res = inner_product(comm_handler, mesh, map, penalisation, is_stag=.false., & u_in=u_inner, & use_abs=.true., & use_vol_avg=.true., & use_full_domain=.true.) fine = almost_equal(res, ref, rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! Test lineout computation function x_coords = [0.26_GP, 0.28_GP, 0.30_GP, 0.32_GP, 0.34_GP] y_coords = [0.0_GP, 0.0_GP, 0.0_GP, 0.0_GP, 0.0_GP] n_points = 5 int_order = 3 reference_lineout(:,1) = [ 0.279273035109860E+000_GP, 0.328639290778531E+000_GP, & 8.273207754419415E-002_GP, -0.237405069737295E+000_GP, & -0.344240215573887E+000_GP] reference_lineout(:,2) = [ 0.558546070219721E+000_GP, 0.657278581557062E+000_GP, & 0.165464155088388E+000_GP, -0.474810139474591E+000_GP, & -0.688480431147774E+000_GP ] reference_lineout(:,3) = [ 0.279273035109861E+000_GP, 0.328639290778531E+000_GP, & 8.273207754419415E-002_GP, -0.237405069737295E+000_GP, & -0.344240215573887E+000_GP ] reference_lineout(:,4) = [ 0_GP, 0_GP, 0_GP, 0_GP, 0_GP ] call interpolate_lineout(mesh, n_points, int_order, x_coords, y_coords, & u_in=u, u_out=res_lineout) do l = 1, 5 fine = almost_equal(res_lineout(l), reference_lineout(l, rank+1), & rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) enddo if (rank == 0) then write(*,'(A80)') & 'test_diagnostics_helpers complete -----------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine end module