init_diffusion_m.f90 Source File


Contents

Source Code


Source Code

module init_diffusion_m
    !! Sets initial state for diffusion model
    use MPI
    use comm_handler_m, only : comm_handler_t
    use parallelisation_setup_m, only : is_rank_info_writer
    use netcdf
    use precision_grillix_m, only : GP
    use parallel_map_m, only : parallel_map_t
    use variable_m, only: variable_t
    use snapshot_m, only : snapshot_t
    use mms_diffusion_m, only : mms_sol_diffusion
    
    ! From PARALLAX
    use error_handling_grillix_m, only: handle_error, error_info_t
    use status_codes_grillix_m, only : GRILLIX_ERR_OTHER, GRILLIX_ERR_NAMELIST
    use screen_io_m, only :  get_stdout
    use equilibrium_m, only : equilibrium_t
    use mesh_cart_m, only : mesh_cart_t  
    implicit none    

    public :: build_initial_state

contains

    subroutine build_initial_state(comm_handler, equi, mesh, map, init_select, filename, snaps, isnaps, tau, var)
        !! Builds initial state
        type(comm_handler_t), intent(in) :: comm_handler
        !! Communicators
        class(equilibrium_t), intent(inout) :: equi
        !! Equilibrium
        type(mesh_cart_t), intent(inout) :: mesh
        !! Mesh within poloidal plane
        type(parallel_map_t), intent(in) :: map
        !! Parallel map
        character(len=*), intent(in) :: init_select
        !! Type of initial state (BLOB, CONTINUE, MMS)
        character(len=*), intent(in) :: filename
        !! Filename, where parameters are possibly read from
        type(snapshot_t), intent(in) :: snaps
        !! Snapshot
        integer, intent(out) :: isnaps
        !! Snapshot number
        real(GP), intent(out) :: tau
        !! Time of initial state
        type(variable_t), intent(inout), target :: var
        !! Variable with initial state

        integer :: io_error
        character(len=256) :: io_errmsg 
        
        real(GP) :: x0, y0, phi0, wx, wy, wphi, amp, wpar
        ! parameters for blob and aligned blob as initial state  

        type(variable_t), allocatable, dimension(:) :: tmp_var
        
        integer :: l
    
        real(GP) :: x, y, phi

        namelist / init_blob / x0, y0, phi0, wx, wy, wphi, amp
        namelist / init_blob_aligned / x0, y0, phi0, wx, wy, wpar, amp


        select case(init_select) 

            case('CONTINUE')
            
                allocate(tmp_var(1))
                call snaps%read_snaps_last(comm_handler, isnaps, tau, tmp_var)
                var = tmp_var(1)
                deallocate(tmp_var)
                                          
                if (is_rank_info_writer) then
                    write(get_stdout(),*)'Continuing simulation at ------------------------------'
                    write(get_stdout(),201)isnaps, tau
201                 FORMAT(1X,'isnaps       = ',I8   /, &
                           1X,'tau          = ',ES14.6E3)  
                    write(get_stdout(),*)'-------------------------------------------------------'  
                endif

            case('BLOB')
                ! blob as initial state
                tau = 0.0_GP
                isnaps = 1

                open(unit = 20, file = filename, 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

                ! some default values                    
                x0 = 0.3_GP
                y0 = 0.0_GP 
                phi0 = 0.0_GP
                wx = 0.03_GP
                wy = 0.03_GP
                wphi = 0.0_GP
                amp = 1.0_GP*map%get_nplanes()
    
                read(20, nml = init_blob, iostat = io_error, iomsg = io_errmsg)
                if (io_error /= 0) then
                    call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
                endif    

                if (is_rank_info_writer) then
                    write(get_stdout(),*)'Blob as initial state ---------------------------------'
                    write(get_stdout(), nml=init_blob)
                    write(get_stdout(),*)'-------------------------------------------------------'  
                endif

                call var%init(comm_handler, ndim = mesh%get_n_points(), staggered = .false.)
                call var%set_blob_toroidal(mesh, x0, y0, phi0, wx, wy, wphi, amp)
                
                close(20)

            case('BLOB_ALIGNED')
                ! Aligned blob as initial state
                tau = 0.0_GP
                isnaps = 1

                open(unit = 20, file = filename, 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

                ! some default values                    
                x0 = 0.3_GP
                y0 = 0.0_GP
                phi0 = 0.0_GP
                wx = 0.03_GP
                wy = 0.03_GP
                wpar = 1.0_GP
                amp = 1.0_GP
    
                read(20, nml = init_blob_aligned, iostat = io_error, iomsg = io_errmsg)
                if (io_error /= 0) then
                    call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
                endif    

                if (is_rank_info_writer) then
                    write(get_stdout(),*)'Aligned blob as initial state -------------------------'
                    write(get_stdout(), nml=init_blob_aligned)
                    write(get_stdout(),*)'-------------------------------------------------------'  
                endif

                call var%init(comm_handler, ndim = mesh%get_n_points(), staggered = .false.)
                call var% set_blob_aligned(equi, mesh, x0, y0, phi0, wx, wy, wpar, amp)
        
                close(20)

             case('MMS')            
                ! Set initial state to MMS solution
                tau = 0.0_GP
                isnaps = 1
                    
                phi = map%get_phi_cano()

                call var%init(comm_handler, mesh%get_n_points(), staggered = .false.)
                !$OMP PARALLEL PRIVATE(l, x, y)
                !$OMP DO
                do l = 1, mesh%get_n_points()
                    x = mesh%get_x(l)
                    y = mesh%get_y(l)
                    var%vals(l) = mms_sol_diffusion(equi, x, y, phi, tau)
                enddo
                !$OMP END DO
                !$OMP END PARALLEL
                
                if (is_rank_info_writer) then
                    write(get_stdout(),*)'MMS as initial state ----------------------------------'
                    write(get_stdout(),*)'-------------------------------------------------------'  
                endif

             case default
                call handle_error('Init_select not valid or not specified', &
                                  GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                                  error_info_t(init_select))
        end select
        
        close(20)
         
    end subroutine

end module