static_data_m.f90 Source File


Contents

Source Code


Source Code

module static_data_m
    !! Holds static data that is time independent.
    use NETCDF
    use technical_constants_m, only : PATHLEN_MAX
    use constants_m, only : TWO_PI
    use screen_io_m, only : get_stdout
    use precision_grillix_m, only : GP, GP_NAN
    use parallelisation_setup_m, only : is_rank_info_writer
    use error_handling_grillix_m, only: handle_error, handle_error_netcdf, error_info_t
    use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST, GRILLIX_ERR_STATIC_DATA, &
        GRILLIX_ERR_CMD, GRILLIX_WRN_GENERAL
    use comm_handler_m, only : comm_handler_t
    use equilibrium_m, only : equilibrium_t
    use mesh_cart_m, only : mesh_cart_t
    use multigrid_m, only: multigrid_t
    use equilibrium_storage_m, only : equilibrium_storage_t
    use helmholtz_solver_m, only : helmholtz_solver_t
    use parallel_map_m, only : parallel_map_t
    use penalisation_m, only : penalisation_t
    use parbnd_taylor_m, only : parbnd_taylor_t
    use polars_m, only : polars_t
    use parallel_target_flux_m, only : parallel_target_flux_t
    use perp_bnd_flux_m, only : perp_bnd_flux_t
    use mask_data_m, only : mask_data_t
    use perf_m, only : perf_start, perf_stop, perf_print
    implicit none

    logical, protected :: static_data_is_initialized
    !! True if full static_data set is initialized

    class(equilibrium_t), public, allocatable :: equi
    !! Equilibrium
    logical, protected :: equi_is_initialized = .false.
    !! True if equi is initialized

    type(mesh_cart_t), public :: mesh_cano
    !! Mesh (canonical) within poloidal plane
    type(mesh_cart_t), public :: mesh_stag
    !! Mesh (staggered) within poloidal plane
    type(multigrid_t), public :: multigrid_cano
    !! Multigrid (canonical)
    type(multigrid_t), public :: multigrid_stag
    !! Multigrid (staggered)
    logical, protected :: multigrid_is_initialized = .false.
    !! True if multigrids and meshes are initialized

    type(equilibrium_storage_t), public, target :: equi_on_cano
    !! Equilibrim quantities on canonical mesh
    type(equilibrium_storage_t), public, pointer :: equi_on_stag => null()
    !! Equilibrim quantities on staggered mesh.
    !! Points to equi_on_cano for axisymmetric case
    !! and to equi_on_stag_mem for non-axisymmetric case
    type(equilibrium_storage_t), private, allocatable, target :: equi_on_stag_mem
    !! Equilibrim quantities on staggered mesh,
    !! only allocated for non-axisymetric case
    logical, protected :: equi_on_mesh_is_initialized = .false.
    !! True, if equi_on_cano/stag are initialized

    class(helmholtz_solver_t), public, allocatable, target :: field_solver_cano
    !! Elliptic 2D field solver on canonical mesh
    class(helmholtz_solver_t), public, pointer :: field_solver_stag => null()
    !! Elliptic 2D field solver on staggered mesh
    !! Points to field_solver_cano for axisymmetric case
    !! and to field_solver_stag_mem for non-axisymmetric case
    class(helmholtz_solver_t), private, allocatable, target :: field_solver_stag_mem
    !! Elliptic 2D solver on staggered mesh
    !! only allocated for non-axisymetric case
    logical, protected :: field_solver_is_initialized = .false.
    !! True, if field_solver is initialized

    type(parallel_map_t), public :: map
    !! Parallel map
    logical, protected :: map_is_initialized = .false.
    !! True, if map is initialized

    type(penalisation_t), public, target :: parbnd_immersed_cano
    !! Immersed boundaries on canonical mesh
    type(penalisation_t), public, pointer :: parbnd_immersed_stag => null()
    !! Immersed boundaries on staggered mesh.
    !! Points to parbnd_immersed_cano for axisymmetric case
    !! and to parbnd_immersed_stag_mem for non-axisymmetric case
    type(penalisation_t), private, allocatable, target :: parbnd_immersed_stag_mem
    !! Immersed boundaries for staggered mesh
    !! only allocated for non-axisymetric case
    logical, protected :: parbnd_immersed_is_initialized = .false.
    !! True, if parbnd_immersed is initialized

    type(parbnd_taylor_t), public, target :: parbnd_taylor_cano
    !! Parallel boundary information based on taylor method for canonical mesh
    type(parbnd_taylor_t), public, pointer :: parbnd_taylor_stag => null()
    !! parbnd_taylor on staggered mesh.
    !! Points to parbnd_taylor_cano for axisymmetric case
    !! and to parbnd_taylor_stag_mem for non-axisymmetric case
    type(parbnd_taylor_t), private, allocatable, target :: parbnd_taylor_stag_mem
    !! parbnd_taylor for staggered mesh
    !! only allocated for non-axisymetric case
    logical, protected :: parbnd_taylor_is_initialized = .false.
    !! True, if parbnd_taylor is initialized

    type(polars_t), public, target :: polars_cano
    !! Polar grid, map and `partial` (flux)-surface operators for canonical mesh
    type(polars_t), public, pointer :: polars_stag => null()
    !! Polar grid, map and `partial` (flux)-surface operators for staggered mesh
    !! Points to polars_cano for axisymmetric case
    !! and to polars_stag_mem for non-axisymmetric case
    type(polars_t), private, allocatable, target:: polars_stag_mem
    !! Polar grid, map and `partial` (flux)-surface operators for staggered mesh
    !! only allocated for non-axisymetric case
    logical, protected :: polars_is_initialized = .false.
    !! True, if polars is initialized

    type(parallel_target_flux_t), public, target :: parflux_utils_cano
    !! Parallel target flux markers and operators for canonical mesh
    type(parallel_target_flux_t), public, pointer :: parflux_utils_stag => null()
    !! Parallel target flux markers and operators for staggered mesh
    !! Points to parflux_utils_cano for axisymmetric case
    !! and to parflux_utils_stag_mem for non-axisymmetric case
    type(parallel_target_flux_t), private, allocatable, target:: parflux_utils_stag_mem
    !! Parallel target flux markers and operators for staggered mesh
    !! only allocated for non-axisymetric case
    type(perp_bnd_flux_t), public, target :: perp_bnd_flux_cano
    !! Perpendicular flux utilities for canonical mesh
    type(perp_bnd_flux_t), public, pointer :: perp_bnd_flux_stag => null()
    !! Perpendicular flux utilities for staggered mesh
    !! Points to perp_bnd_flux_cano for axisymmetric case
    !! and to perp_bnd_flux_stag_mem for non-axisymmetric case
    type(perp_bnd_flux_t), private, allocatable, target:: perp_bnd_flux_stag_mem
    !! Perpendicular flux utilities for staggered mesh
    !! only allocated for non-axisymetric case
    logical, protected :: bnd_flux_is_initialized = .false.
    !! True, if bnd_flux utilities are initialized

    type(mask_data_t), public, target :: masks_cano
    !! Mask arrays for canonical mesh
    type(mask_data_t), public, pointer :: masks_stag => null()
    !! Mask arrays for staggered mesh
    !! Points to masks_cano for axisymmetric case
    !! and to masks_stag_mem for non-axisymmetric case
    type(mask_data_t), private, allocatable, target:: masks_stag_mem
    !! Mask arrays for staggered mesh
    !! only allocated for non-axisymetric case
    logical, protected :: masks_is_initialized = .false.
    !! True, if masks are initialized

    enum, bind(c)
        !! Available values to control build and NetCDF I/O of individual components
        enumerator :: NONE = 0
        !! Does not build component
        enumerator :: BUILD_AND_WRITE = 1
        !! Builds component and writes it to NetCDF file
        enumerator :: READ = 2
        !! Reads component from netdcf file
        enumerator :: BUILD_NOWRITE = 3
        !! Builds component without NetCDF output
    end enum

    ! Parameters that control build workflow
    integer, protected :: cntrl_build_multigrid = BUILD_AND_WRITE
    !! Controls build and NetCDF I/O for multigrid/mesh
    integer, protected :: cntrl_build_equi_on_mesh = BUILD_AND_WRITE
    !! Controls build and NetCDF I/O for equilibrium storage on mesh
    integer, protected :: cntrl_build_map = BUILD_AND_WRITE
    !! Controls build and NetCDF I/O for map
    integer, protected :: cntrl_build_parbnd_immersed = BUILD_AND_WRITE
    !! Controls build and NetCDF I/O for parbnd_immersed
    integer, protected :: cntrl_build_parbnd_taylor = NONE
    !! Controls build and NetCDF I/O for parbnd_taylor
    integer, protected :: cntrl_build_polars = BUILD_AND_WRITE
    !! Controls build and NetCDF I/O for polars
    integer, protected :: cntrl_build_bnd_flux = NONE
    !! Controls build and NetCDF I/O for bnd_flux utilities
    integer, protected :: cntrl_build_masks = BUILD_AND_WRITE
    !! Controls build and NetCDF I/O for masks
    integer, protected :: cntrl_build_vtkmesh = NONE
    !! Controls build for vtk (3D visualisation) mesh
    !! Only available, if compiled with -DPARALLAX_ENABLE_VTK=ON
    character(len=PATHLEN_MAX), protected :: static_data_dirpath = 'trunk'
    !! Directory path where to read and/or write NetCDF files from
    character(len=PATHLEN_MAX), protected :: static_data_parpath = 'params_static_data.nml'
    !! Filepath where to read parameters for initialisation of static data
    logical, protected :: enable_debug_output = .true.
    !! Set to true, to activate more verbose output

    public :: static_all_components_init
    public :: finalize_static
    public :: read_build_static
    public :: display_build_static

    public :: static_equi_init
    public :: static_multigrid_init
    public :: static_equi_on_mesh_init
    public :: static_field_solver_init
    public :: static_map_init
    public :: static_parbnd_immersed_init
    public :: static_polars_init
    public :: static_bnd_flux_init

    interface

        module subroutine static_equi_init(par_filepath)
            !! Initializes equi as part of static data
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read equilibrium parameters from
        end subroutine

        module subroutine static_multigrid_init(comm_handler, par_filepath, nf90_filepath)
            !! Initializes multigrid and mesh (canonical and staggered) as part of static data
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read multigrid parameters from
            character(len=*), intent(in), optional :: nf90_filepath
            !! Filename, where to read/write multigrid NetCDF file,
            !! If not present NetCDF I/O happens to static_data_dirpath directory
        end subroutine

        module subroutine static_equi_on_mesh_init(comm_handler, nf90_filepath)
            !! Initializes equi_on_mesh (canonical and staggered) as part of static data
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: nf90_filepath
            !! Filename, where to read/write equi_on_mesh NetCDF file,
            !! If not present NetCDF I/O happens to static_data_dirpath directory
        end subroutine

        module subroutine static_field_solver_init(par_filepath)
            !! Initializes field_solver (canonical and staggered) as part of static data
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read field_solver parameters from
        end subroutine

        module subroutine static_map_init(comm_handler, par_filepath, nf90_filepath)
            !! Initializes parallel map as part of static data
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read map parameters from
            character(len=*), intent(in), optional :: nf90_filepath
            !! Filename, where to read/write map NetCDF file,
            !! If not present NetCDF I/O happens to static_data_dirpath directory
        end subroutine

        module subroutine static_parbnd_immersed_init(comm_handler, par_filepath, nf90_filepath)
            !! Initializes parbnd_immersed as part of static data
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read parbnd_immersed parameters from
            character(len=*), intent(in), optional :: nf90_filepath
            !! Filename, where to read/write parbnd_immersed NetCDF file,
            !! If not present NetCDF I/O happens to static_data_dirpath directory
        end subroutine

        module subroutine static_parbnd_taylor_init(comm_handler, par_filepath, nf90_filepath)
            !! Initializes parbnd_taylor as part of static data
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read parbnd_taylor parameters from
            character(len=*), intent(in), optional :: nf90_filepath
            !! Filename, where to read/write parbnd_taylor NetCDF file,
            !! If not present NetCDF I/O happens to static_data_dirpath directory
        end subroutine

        module subroutine static_polars_init(comm_handler, par_filepath, nf90_filepath)
            !! Initializes polars as part of static data
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read polars parameters from
            character(len=*), intent(in), optional :: nf90_filepath
            !! Filename, where to read/write polars NetCDF file,
            !! If not present NetCDF I/O happens to static_data_dirpath directory
        end subroutine

        module subroutine static_bnd_flux_init(comm_handler, par_filepath, nf90_filepath)
            !! Initializes bnd_flux as part of static data
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read bnd_flux parameters from
            character(len=*), intent(in), optional :: nf90_filepath
            !! Filename, where to read/write bnd_flux NetCDF file,
            !! If not present NetCDF I/O happens to static_data_dirpath directory
        end subroutine

        module subroutine static_masks_init(comm_handler, nf90_filepath)
            !! Initializes masks as part of static data
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: nf90_filepath
            !! Filename, where to read/write masks NetCDF file,
            !! If not present NetCDF I/O happens to static_data_dirpath directory
        end subroutine

        module subroutine build_vtk_mesh(comm_handler, par_filepath, dirpath)
            !! Builds VTK mesh for 3D visualisation
            !! Functionality only available if compiled with -DENABLE_VTK=ON
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in), optional :: par_filepath
            !! Filename, where to read VTK parameters from
            character(len=*), intent(in), optional :: dirpath
            !! Directory, where to write VTK meshes to
        end subroutine

    end interface

contains

    subroutine static_all_components_init(comm_handler, par_filepath)
        !! Initializes a complete static data set with all components
        type(comm_handler_t), intent(in) :: comm_handler
        !! Communication handler
        character(len=*), intent(in), optional :: par_filepath
        !! Filename, where to read parameters for all components from
        integer :: ierr, cmd_exit, cmd_err
        logical :: any_netcdf_output_requested, dir_exists

        ! Read control parameters for build of static data
        if (present(par_filepath)) then
            call read_build_static(par_filepath)
        else
            call read_build_static()
        endif
        call display_build_static()

        ! Check if any NetCDF output of static data is requested
        if ( (cntrl_build_multigrid == BUILD_AND_WRITE) .or. &
             (cntrl_build_equi_on_mesh == BUILD_AND_WRITE) .or. &
             (cntrl_build_map == BUILD_AND_WRITE) .or. &
             (cntrl_build_parbnd_immersed == BUILD_AND_WRITE) .or. &
             (cntrl_build_polars == BUILD_AND_WRITE) .or. &
             (cntrl_build_bnd_flux == BUILD_AND_WRITE)) then
            any_netcdf_output_requested = .true.
        endif

        ! Create directory, where NetCDF data is is written
        if ((any_netcdf_output_requested) .and. (is_rank_info_writer)) then
            call execute_command_line('mkdir -p ' // trim(static_data_dirpath), &
                                        exitstat=cmd_exit, cmdstat=cmd_err)
            if (cmd_err /= 0) then
                call handle_error("Command execution failed: mkdir -p " // trim(static_data_dirpath), &
                                    GRILLIX_ERR_CMD, __LINE__, __FILE__, &
                                    error_info_t('cmd_err: ',[cmd_err]))
            endif
            if (cmd_exit /= 0) then
                call handle_error("Command execution returned with warning: mkdir -p" &
                                    // trim(static_data_dirpath), &
                                    GRILLIX_ERR_CMD, __LINE__, __FILE__, &
                                    error_info_t('cmd_exit: ',[cmd_exit]))
            endif
        endif
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        ! Initialize components in consistent order,
        ! the MPI barriers ensure that NetCDF output is complete before next component is built
        call perf_start("static_data/equi_init")
        call static_equi_init(par_filepath=static_data_parpath)
        call perf_stop("static_data/equi_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/multigrid_init")
        call static_multigrid_init(comm_handler, par_filepath=static_data_parpath)
        call perf_stop("static_data/multigrid_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/equi_on_mesh_init")
        call static_equi_on_mesh_init(comm_handler)
        call perf_stop("static_data/equi_on_mesh_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/field_solver_init")
        call static_field_solver_init(par_filepath=static_data_parpath)
        call perf_stop("static_data/field_solver_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/map_init")
        call static_map_init(comm_handler, par_filepath=static_data_parpath)
        call perf_stop("static_data/map_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/parbnd_immersed_init")
        call static_parbnd_immersed_init(comm_handler, par_filepath=static_data_parpath)
        call perf_stop("static_data/parbnd_immersed_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/parbnd_taylor_init")
        call static_parbnd_taylor_init(comm_handler, par_filepath=static_data_parpath)
        call perf_stop("static_data/parbnd_taylor_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/polars_init")
        call static_polars_init(comm_handler, par_filepath=static_data_parpath)
        call perf_stop("static_data/polars_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/bnd_flux_init")
        call static_bnd_flux_init(comm_handler, par_filepath=static_data_parpath)
        call perf_stop("static_data/bnd_flux_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/masks_init")
        call static_masks_init(comm_handler)
        call perf_stop("static_data/masks_init")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        call perf_start("static_data/build_vtk_mesh")
        call build_vtk_mesh(comm_handler, par_filepath=static_data_parpath)
        call perf_stop("static_data/build_vtk_mesh")
        call MPI_barrier(comm_handler%get_comm_cart(), ierr)

        ! Wrap up
        call finalize_static()

        if (enable_debug_output) then
            ! Performance statics summed over all ranks
            call perf_print(comm=comm_handler%get_comm_cart())
        endif

    end subroutine

    subroutine finalize_static()
        !! Cleans up temporary data, that was used during initialization stage,
        !! and sets static_data_is_initialized to true

        call mesh_cano%deallocate_shaped_cart_arrays()
        call mesh_stag%deallocate_shaped_cart_arrays()

        static_data_is_initialized = .true.

    end subroutine

    subroutine read_build_static(par_filepath)
        !! Reads parameters that control build of static data
        character(len=*), intent(in), optional :: par_filepath
        !! Filename, where to read control build parameters from
        character(len=32) :: build_multigrid = 'BUILD_AND_WRITE'
        character(len=32) :: build_equi_on_mesh = 'BUILD_AND_WRITE'
        character(len=32) :: build_map = 'BUILD_AND_WRITE'
        character(len=32) :: build_parbnd_immersed = 'BUILD_AND_WRITE'
        character(len=32) :: build_parbnd_taylor = 'NONE'
        character(len=32) :: build_polars = 'BUILD_AND_WRITE'
        character(len=32) :: build_bnd_flux = 'NONE'
        character(len=32) :: build_masks = 'BUILD_AND_WRITE'
        character(len=32) :: build_vtkmesh = 'NONE'

        namelist / build_static / build_multigrid, build_equi_on_mesh, build_map, &
            build_parbnd_immersed, build_parbnd_taylor, build_polars, build_bnd_flux, &
            build_masks, build_vtkmesh, static_data_dirpath, enable_debug_output

        integer :: io_unit, io_error
        character(len=256) :: io_errmsg

        if (present(par_filepath)) then
            static_data_parpath = par_filepath
        endif

        ! Read build_static parameters
        open(newunit=io_unit, file=static_data_parpath, status='old', action='read', iostat=io_error, iomsg=io_errmsg)
        if (io_error /= 0) then
            call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
        endif
        read(io_unit, nml=build_static, iostat=io_error, iomsg=io_errmsg)
        if (io_error /= 0) then
            call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
        endif
        close(io_unit)

        cntrl_build_multigrid = convert_build_string_to_enum(build_multigrid)
        cntrl_build_equi_on_mesh = convert_build_string_to_enum(build_equi_on_mesh)
        cntrl_build_map = convert_build_string_to_enum(build_map)
        cntrl_build_parbnd_immersed = convert_build_string_to_enum(build_parbnd_immersed)
        cntrl_build_parbnd_taylor = convert_build_string_to_enum(build_parbnd_taylor)
        cntrl_build_polars = convert_build_string_to_enum(build_polars)
        cntrl_build_bnd_flux = convert_build_string_to_enum(build_bnd_flux)
        cntrl_build_masks = convert_build_string_to_enum(build_masks)
        cntrl_build_vtkmesh = convert_build_string_to_enum(build_vtkmesh)

    contains
        function convert_build_string_to_enum(build_string) result(res)
            !! Converts build string to an emumerator, as defined in the header above
            character(len=*), intent(in) :: build_string
            integer :: res
            select case(trim(adjustl(build_string)))
            case("BUILD_AND_WRITE")
                res = BUILD_AND_WRITE
            case("BUILD_NOWRITE")
                res = BUILD_NOWRITE
            case("READ")
                res = READ
            case("NONE")
                res = NONE
            case default
                call handle_error('build_string not valid:' // build_string, &
                                   GRILLIX_ERR_STATIC_DATA, __LINE__, __FILE__)
            end select
        end function
    end subroutine

    subroutine display_build_static()
        !! Writes control build parameters of static data parameters to screen

        namelist / cntrl_build_static / cntrl_build_multigrid, cntrl_build_equi_on_mesh, cntrl_build_map, &
            cntrl_build_parbnd_immersed, cntrl_build_parbnd_taylor, cntrl_build_polars, cntrl_build_bnd_flux, &
            cntrl_build_masks, cntrl_build_vtkmesh, static_data_parpath, static_data_dirpath, enable_debug_output

        if (is_rank_info_writer) then
            write(get_stdout(), nml=cntrl_build_static)
        endif
    end subroutine

end module