commons_misc_m.f90 Source File


Contents

Source Code


Source Code

module commons_misc_m
    !! Miscellaneous routines that are common to all models
    use MPI
    use comm_handler_m, only : comm_handler_t
    use screen_io_m, only :  get_stdout
    use precision_grillix_m, only : GP, MPI_GP

    public :: error_ananum
    
contains

    subroutine error_ananum(comm_handler, ndim, ana, num, nrm2, nrmsup, err2, errsup)
        !! Computes the norm and difference between two fields, useful for MMS-analysis
        type(comm_handler_t), intent(in) :: comm_handler
        !! Communicators
        integer, intent(in) :: ndim
        !! Number of points
        real(GP), dimension(ndim), intent(in) :: ana
        !! Field 1, e.g. analytic MMS-solution
        real(GP), dimension(ndim), intent(in) :: num
        !! Field 2, e.g. numerical solution
        real(GP), intent(out) :: nrm2
        !! L2 norm of Field 1
        real(GP), intent(out) :: nrmsup
        !! Supremums norm of Field 1
        real(GP), intent(out) :: err2
        !! Error in L2 norm
        real(GP), intent(out) :: errsup
        !! Error in Supremumsnorm

        integer :: i, ierr
        real(GP) :: diff, nrm_points

        nrm2     = 0.0_GP
        nrmsup   = 0.0_GP 
        err2     = 0.0_GP
        errsup   = 0.0_GP

        !$OMP PARALLEL PRIVATE(i, diff)
        !$OMP DO REDUCTION(+: nrm2, err2) Reduction(max: nrmsup, errsup)
        do i = 1, ndim

            diff = ana(i) - num(i)
  
            nrm2   = nrm2 + ana(i)**2
            nrmsup = max(nrmsup, abs(ana(i)))
            err2   = err2 + diff**2
            errsup = max(errsup, abs(diff))
        enddo
        !$OMP END DO
        !$OMP END PARALLEL 

        call MPI_allreduce(MPI_IN_PLACE, nrm2,   1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr)
        call MPI_allreduce(MPI_IN_PLACE, nrmsup, 1, MPI_GP, MPI_MAX, comm_handler%get_comm_planes(), ierr)
        call MPI_allreduce(MPI_IN_PLACE, err2,   1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr)
        call MPI_allreduce(MPI_IN_PLACE, errsup, 1, MPI_GP, MPI_MAX, comm_handler%get_comm_planes(), ierr)
        
        nrm_points = 1.0_GP*ndim*comm_handler%get_nplanes()

        nrm2 = sqrt( nrm2 / nrm_points )
        err2 = sqrt( err2 / nrm_points )     
            
    end subroutine

end module