test_apar_fluct_evaluate_m.pf Source File

module test_apar_fluct_evaluate_m
    !! Tests regarding apar_shift module
    !! Test evaluation of apar_fluct in an axisymmetric geomtery
    !! which involves toroidal average
    use pfunit
    use MPI
    use funit
    use comm_handler_m, only : comm_handler_t
    use precision_grillix_m, only : GP
    use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST
    use error_handling_grillix_m, only: handle_error
    use helpers_m, only : almost_equal
    use helpers_mesh_sample_m, only : mesh_sample_circular
    use screen_io_m, only : get_stdout
    use model_apar_shift_m, only : apar_shift_t
    use commons_misc_m, only : error_ananum
    use equilibrium_m, only: equilibrium_t
    use mesh_cart_m, only : mesh_cart_t
    use variable_m, only : variable_t
    use constants_m, only : two_pi, pi
    implicit none
contains

    @test(npes = [4])
    subroutine test_apar_fluct_evaluate(this)
        !! Tests evaluate_apar_fluct function 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_stag
        type(apar_shift_t) :: ashift
        type(variable_t) :: apar
        type(variable_t) :: apar_fluct
        type(variable_t) :: apar_fluct_ref
        
        real(GP) :: phi_stag, phi_cano, dphi
        logical :: fine
        real(GP) :: nrm2, nrmsup, err2, errsup
        integer :: l

        comm_world = this%getMpiCommunicator()
        call MPI_comm_rank(comm_world, rank, ierr)

        if (rank == 0) then
            write(*,*)''
            write(*,*)''
            write(*,'(A80)')'test_apar_fluct_evaluate ---------------------------------------------------'
        endif

        call comm_handler%init(comm_world, 4, 1)
        
        ! Create samplemesh ---------------------------------------------------------------------------
        call comm_handler%init(comm_world, 4, 1)
        dphi = TWO_PI / comm_handler%get_nplanes() ! Pi/2
        phi_cano = dphi * comm_handler%get_plane() ! 0, Pi/2, Pi, 3Pi/2
        phi_stag = phi_cano + dphi / 2.0_GP ! Pi/4, 3Pi/4, 5Pi/4, 7Pi/4

        ! Create samplemeshes -------------------------------------------------------------------------
        call mesh_sample_circular(equi, phi_stag, mesh_stag)

        ! Initialising ---------------------------------------------------------------------
        call apar%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.)
        call ashift%apar_shift%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.)       
        call apar_fluct%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.)    
        call apar_fluct_ref%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.)    

        !$omp parallel default(none) private(l) &
        !$omp shared(mesh_stag, apar, ashift, apar_fluct_ref, phi_stag) 
        !$omp do
        do l = 1, mesh_stag%get_n_points()
            apar%vals(l) = cos(phi_stag) + 1.0_GP
            ashift%apar_shift%vals(l) = 9.0_GP * sin(phi_stag) + 4.0_GP * sin(phi_stag)**2 + phi_stag / pi
            apar_fluct_ref%vals(l) = cos(phi_stag) - 2.0_GP
        enddo
        !$omp end do
        !$omp end parallel
        call ashift%eval_fluct(comm_handler, equi, mesh_stag, apar, apar_fluct)
      
        ! Check results -------------------------------------------------------        
        call error_ananum(comm_handler, mesh_stag%get_n_points(), &
                          apar_fluct_ref%vals, apar_fluct%vals, &
                          nrm2, nrmsup, err2, errsup)
        fine = almost_equal(err2, 0.0_GP,  0.0_GP, 1.0E-9_GP)
        @assertTrue(fine)

        if (rank == 0) then
            write(*,'(A80)')'test_apar_fluct_evaluate complete ----------------------------'
            write(*,*)''
            write(*,*)''
        endif

    end subroutine
    
end module