helmholtz_netcdfio_m.f90

Source Code

module helmholtz_netcdfio_m
    !! I/O to NETCDF format for a Helmholtz problem,
    !! useful for debugging and benchmarks
    use netcdf
    use precision_m, only : FP
    use screen_io_m, only : get_stderr, nf90_handle_err
    use comm_handling_m, only: is_master
    use descriptors_m, only : DISTRICT_CORE, DISTRICT_WALL, &
                              DISTRICT_DOME, DISTRICT_OUT
    use csrmat_m, only : csrmat_t
    use mesh_cart_m, only : mesh_cart_t
    use helmholtz_operator_m, only : build_helmholtz_csr
    implicit none

    public :: read_netcdf_helmholtz
    public :: write_netcdf_helmholtz
contains

    subroutine write_netcdf_helmholtz(fgid, mesh, bnd_type_core, &
                                      bnd_type_wall, bnd_type_dome, &
                                      bnd_type_out, &
                                      co, lambda, xi, &
                                      rhs, guess, sol, &
                                      hcsr_write_on)
        !! Writes a Helmholtz problem to NETCDF file
        integer, intent(in) :: fgid
        !! NETCDF file or group id
        type(mesh_cart_t), intent(in) :: mesh
        !! Mesh
        integer, intent(in) :: bnd_type_core
        !! Boundary descriptor for core boundary
        integer, intent(in) :: bnd_type_wall
        !! Boundary descriptor for wall boundary
        integer, intent(in) :: bnd_type_dome
        !! Boundary descriptor for dome boundary
        integer, intent(in) :: bnd_type_out
        !! Boundary descriptor for outer(mask) boundary
        real(FP), dimension(mesh%get_n_points()), intent(in) :: co
        !! Coefficient within Helmholtz operator
        real(FP), dimension(mesh%get_n_points_inner()), intent(in) :: lambda
        !! Coefficient within Helmholtz operator
        real(FP), dimension(mesh%get_n_points_inner()), intent(in) :: xi
        !! Coefficient within Helmholtz operator
        real(FP), dimension(mesh%get_n_points()), intent(in) :: rhs
        !! Right hand side
        real(FP), dimension(mesh%get_n_points()), intent(in) :: guess
        !! Initial guess
        real(FP), dimension(mesh%get_n_points()), intent(in) :: sol
        !! Solution
        logical, intent(in) :: hcsr_write_on
        !! Switch to write also Helmoltz matrix in CSR format
        !! Useful for debugging, benchmarking wwith external tools

        integer :: iddim_n_points, iddim_n_points_inner, &
                   iddim_n_points_boundary, iddim_n_points_ghost
        integer :: nf90_stat, id_co, id_lambda, id_xi, id_rhs, id_sol, id_guess
        integer :: grpid_hcsr
        integer :: kb, l, district
        integer, dimension(mesh%get_n_points_boundary()) :: bnd_descrs
        real(FP), dimension(mesh%get_n_points()) :: hdiag_inv
        type(csrmat_t) :: hcsr
        character(len=*), parameter :: cname = "write_netcdf_helmholtz"
        ! Name used for logging errors

        ! Write bnd_types as attributes
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'bnd_type_core', &
                                 bnd_type_core)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'bnd_type_wall', &
                                 bnd_type_wall)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'bnd_type_dome', &
                                 bnd_type_dome)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'bnd_type_out', &
                                 bnd_type_out)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Define dimensions
        nf90_stat = nf90_def_dim(fgid, 'n_points', &
                                 mesh%get_n_points(), iddim_n_points)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_dim(fgid, 'n_points_inner', &
                                 mesh%get_n_points_inner(), &
                                 iddim_n_points_inner)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_dim(fgid, 'n_points_boundary', &
                                 mesh%get_n_points_boundary(), &
                                 iddim_n_points_boundary)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_dim(fgid, 'n_points_ghost', &
                                 mesh%get_n_points_ghost(), &
                                 iddim_n_points_ghost)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Define variables
        nf90_stat = nf90_def_var(fgid, 'co', NF90_DOUBLE, &
                                 iddim_n_points, id_co)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'lambda', NF90_DOUBLE, &
                                 iddim_n_points_inner, id_lambda)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'xi', NF90_DOUBLE, &
                                 iddim_n_points_inner, id_xi)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'rhs', NF90_DOUBLE, &
                                 iddim_n_points, id_rhs)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'sol', NF90_DOUBLE, &
                                 iddim_n_points, id_sol)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'guess', NF90_DOUBLE, &
                                 iddim_n_points, id_guess)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Define group for hcsr matrix
        if (hcsr_write_on) then
            nf90_stat = nf90_def_grp(fgid, 'hcsr', grpid_hcsr)
            if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        endif

        ! End of definition
        nf90_stat = nf90_enddef(fgid)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)


        ! Put variables
        nf90_stat = nf90_put_var(fgid, id_co, co)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_lambda, lambda)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_xi, xi)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_rhs, rhs)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_sol, sol)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_guess, guess)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Write hcsr matrix
        if (hcsr_write_on) then

            ! Build matrix
            !$OMP PARALLEL PRIVATE(kb, l, district)
            !$OMP DO
            do kb = 1, mesh%get_n_points_boundary()
                l = mesh%boundary_indices(kb)
                district = mesh%get_district(l)
                if (district == DISTRICT_CORE) then
                    bnd_descrs(kb) = bnd_type_core
                elseif (district == DISTRICT_WALL) then
                    bnd_descrs(kb) = bnd_type_wall
                elseif (district == DISTRICT_DOME) then
                    bnd_descrs(kb) = bnd_type_dome
                elseif (district == DISTRICT_OUT) then
                    bnd_descrs(kb) = bnd_type_out
                else
                    write(get_stderr(),*)'error(write_netcdf_helmholtz): &
                                   district not valid', l, district
                    error stop
                endif
            enddo
            !$OMP END DO
            !$OMP END PARALLEL

            call build_helmholtz_csr(mesh, hcsr, hdiag_inv, &
                                     co, xi, lambda, bnd_descrs)

            call hcsr%write_netcdf(grpid_hcsr)

        endif

    end subroutine

    subroutine read_netcdf_helmholtz(fgid, mesh, bnd_type_core, &
                                     bnd_type_wall, bnd_type_dome, &
                                     bnd_type_out, &
                                     co, lambda, xi, &
                                     rhs, guess, sol)
        !! Reads a Helmholtz problem to NETCDF file
        integer, intent(in) :: fgid
        !! NETCDF file or group id
        type(mesh_cart_t), intent(in) :: mesh
        !! Mesh
        integer, intent(out) :: bnd_type_core
        !! Boundary descriptor for core boundary
        integer, intent(out) :: bnd_type_wall
        !! Boundary descriptor for wall boundary
        integer, intent(out) :: bnd_type_dome
        !! Boundary descriptor for dome boundary
        integer, intent(out) :: bnd_type_out
        !! Boundary descriptor for outer(mask) boundary
        real(FP), dimension(mesh%get_n_points()), intent(out) :: co
        !! Coefficient within Helmholtz operator
        real(FP), dimension(mesh%get_n_points_inner()), intent(out) :: lambda
        !! Coefficient within Helmholtz operator
        real(FP), dimension(mesh%get_n_points_inner()), intent(out) :: xi
        !! Coefficient within Helmholtz operator
        real(FP), dimension(mesh%get_n_points()), intent(out) :: rhs
        !! Right hand side
        real(FP), dimension(mesh%get_n_points()), intent(out) :: guess
        !! Initial guess
        real(FP), dimension(mesh%get_n_points()), intent(out) :: sol
        !! Solution


        integer :: nf90_stat, id_co, id_lambda, id_xi, id_rhs, id_sol, id_guess
        character(len=*), parameter :: cname = "read_netcdf_helmholtz"
        ! Name used for logging errors

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'bnd_type_core', &
                                 bnd_type_core)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'bnd_type_wall', &
                                 bnd_type_wall)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'bnd_type_dome', &
                                 bnd_type_dome)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'bnd_type_out', &
                                 bnd_type_out)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Read co
        nf90_stat = nf90_inq_varid(fgid, 'co', id_co)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_co, co)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Read lambda
        nf90_stat = nf90_inq_varid(fgid, 'lambda', id_lambda)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_lambda, lambda)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Read xi
        nf90_stat = nf90_inq_varid(fgid, 'xi', id_xi)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_xi, xi)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Read rhs
        nf90_stat = nf90_inq_varid(fgid, 'rhs', id_rhs)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_rhs, rhs)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Read sol
        nf90_stat = nf90_inq_varid(fgid, 'sol', id_sol)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_sol, sol)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Read guess
        nf90_stat = nf90_inq_varid(fgid, 'guess', id_guess)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_guess, guess)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

    end subroutine

end module
helmholtz_netcdfio_m.f90 Source File

Contents

Modules

Source Code

Source Code
