mod_gimli.t

Go to the documentation of this file.

!> Module for reading in Legolas .ldat data as initial condition
module mod_gimli
    use, intrinsic :: iso_fortran_env
    use mod_global_parameters
    use mod_hd_phys, only: rho_hd=>rho_, mom_hd=>mom, p_hd=>p_
    use mod_mhd_phys, only: rho_mhd=>rho_, mom_mhd=>mom, p_mhd=>p_ 
    use mod_functions_bfield, only: mag_mhd=>mag
    use mod_eos, only: eos
    use mod_comm_lib, only: mpistop
 
    implicit none
    public
 
    integer, parameter :: dp = real64
    integer :: idirmin0, idirmin
 
    complex(dp), parameter :: ic = (0.0d0, 1.0d0)
 
    integer :: ef_gridpts
    real(dp) :: k2, k3
    integer :: mhd_bool
    real(dp), allocatable :: ef_grid(:)
    complex(dp), allocatable :: rho(:)
    complex(dp), allocatable :: v1(:)
    complex(dp), allocatable :: v2(:)
    complex(dp), allocatable :: v3(:)
    complex(dp), allocatable :: p(:)
    complex(dp), allocatable :: b1(:)
    complex(dp), allocatable :: b2(:)
    complex(dp), allocatable :: b3(:)
 
    integer :: rho_idx, p_idx
    integer, allocatable :: mom_idx(:), mag_idx(:)
 
contains
 
    subroutine read_legolas_data(legolas_file, file_id)
        character(len=100), intent(in) :: legolas_file
        integer, intent(in) :: file_id
 
        open( &
        unit=file_id+mype, &
        file=legolas_file, &
        form='unformatted' &
        )
 
        ! need to read physics_type because mod_physics not yet loaded at this point
        read(file_id+mype) mhd_bool
        read(file_id+mype) ef_gridpts
        read(file_id+mype) k2, k3
 
        call allocate_arrays(ef_gridpts)
        read(file_id+mype) ef_grid
        read(file_id+mype) rho
        read(file_id+mype) v1
        read(file_id+mype) v2
        read(file_id+mype) v3
        read(file_id+mype) p
        if (mhd_bool == 1) then
            read(file_id+mype) b1
            read(file_id+mype) b2
            read(file_id+mype) b3
        end if
 
        if (mhd_bool == 1) then 
            read(file_id+mype) unit_length, unit_numberdensity, unit_temperature, &
            unit_density, unit_pressure, unit_velocity, unit_magneticfield, unit_time
        else 
            read(file_id+mype) unit_length, unit_numberdensity, unit_temperature, &
            unit_density, unit_pressure, unit_velocity, unit_time
        end if
 
        close(file_id+mype)
    end subroutine read_legolas_data
 
    subroutine allocate_arrays(gridpts)
        integer, intent(in) :: gridpts
 
        allocate(ef_grid(gridpts))
        allocate(rho(gridpts))
        allocate(v1, v2, p, mold=rho)
        if (ndim >= 3) allocate(v3, mold=rho)
        if (mhd_bool == 1) then
            allocate(b1, b2, mold=rho)
            if (ndim >= 3) allocate(b3, mold=rho)
        end if
    end subroutine allocate_arrays
 
    subroutine add_perturbation_to_w_array(ixI^L, ixO^L, w, w_index, x)
        integer, intent(in)     :: ixI^L, ixO^L
        real(dp), intent(inout) :: w(ixI^S, nw)
        integer, intent(in)     :: w_index
        real(dp), intent(in)    :: x(ixI^S, ndim)
        complex(dp) :: amplitude, exp_factor(ixI^S), quantity, values(ef_gridpts)
        integer  :: idx^D
        real(dp) :: k1 = 0.0d0
 
        call w_index_to_array(w_index, values)
        exp_factor = {exp(ic * k^d * x(ixi^s, ^d))|*}
 
        {do idx^d = iximin^d, iximax^d|\}
        call ef_amplitude(x(idx^d, 1), ef_grid, values, amplitude)
        quantity = amplitude * exp_factor(idx^d)
        w(idx^d, w_index) = w(idx^d, w_index) + quantity%re
        {end do|\}
    end subroutine add_perturbation_to_w_array
 
    subroutine init_gimli_indices()
        if (mhd_bool == 1) then
            rho_idx = rho_mhd
            p_idx   = p_mhd
 
            if (.not. allocated(mom_idx)) allocate(mom_idx(size(mom_mhd)))
            mom_idx = mom_mhd
 
            if (.not. allocated(mag_idx)) allocate(mag_idx(size(mag_mhd)))
            mag_idx = mag_mhd
        else
            rho_idx = rho_hd
            p_idx   = p_hd
 
            if (.not. allocated(mom_idx)) allocate(mom_idx(size(mom_hd)))
            mom_idx = mom_hd
 
        end if
    end subroutine init_gimli_indices
 
    subroutine w_index_to_array(w_index, array)
        integer, intent(in) :: w_index
        complex(dp), intent(inout) :: array(ef_gridpts)
 
        call init_gimli_indices()
 
        if (w_index == rho_idx) then
            array = rho
        else if (w_index == mom_idx(1)) then
            array = v1
        else if (w_index == mom_idx(2)) then
            array = v2
        end if
 
        if (ndim >= 3) then        
            if (w_index == mom_idx(3)) then
                array = v3
            end if
        end if
        if (w_index == p_idx) then
            array = p
        end if
        if (mhd_bool == 1) then
            if (w_index == mag_idx(1)) then
                array = b1
            else if (w_index == mag_idx(2)) then
                array = b2
            end if
            if (ndim >= 3) then
                if (w_index == mag_idx(3)) then
                    array = b3
                end if
            end if
        end if 
    end subroutine w_index_to_array
 
    subroutine ef_amplitude(x, grid, array, amplitude)
        use mod_geometry, only: coordinate, cylindrical
        real(dp), intent(in) :: x, grid(ef_gridpts)
        complex(dp), intent(in) :: array(ef_gridpts)
        complex(dp), intent(out) :: amplitude
        integer :: idl, idu
        real(dp) :: alpha
 
        if (x <= grid(1) .and. xprobmin1 <= x) then
            amplitude = array(1)
            
            if (coordinate==cylindrical .and. grid(1) > 0) then
                if (k2 .ne. 0.0d0) then
                    alpha = abs(k2)
                else
                    alpha = 2.0d0
                end if
                amplitude = array(1) * (x / grid(1))**alpha
            end if
        else if (x >= grid(size(grid)) .and. x <= xprobmax1) then
            amplitude = array(size(grid))
        else if (x < xprobmin1 .or. x > xprobmax1) then
            amplitude = 0.d0
        else
            idl = maxloc(grid, mask=(grid < x), dim=1)
            idu = minloc(grid, mask=(grid > x), dim=1)
            amplitude = array(idl) + (x - grid(idl)) * &
                (array(idu) - array(idl)) / (grid(idu) - grid(idl))
        end if
    end subroutine ef_amplitude
 
    ! subroutine custom analytics log file
    ! called in "usr_print_log => analytics_log" in mod_usr.t
    subroutine analytics_log
        use mod_global_parameters
        use mod_input_output, only: get_volume_average_func, printlog_default
 
        integer, parameter :: n_modes = 2
        integer, parameter :: my_unit = 123
        character(len=80)  :: fmt_string = '(9(es12.4))', filename
        logical, save      :: alive, visited=.false.
        double precision   :: volume, magn_avg
        double precision   :: Tmax, Tmin, vmax, B1max, B2max, B3max
 
        ! First output the standard log file:
        call printlog_default
 
        ! Now make the custom _c.log file:
        call get_minmax_temperature(tmax,tmin)
        call get_max_velocity(vmax)
        if (mhd_bool == 1) then
            call get_max_b(b1max, b2max, b3max)
            call get_volume_average_func(magnetic, magn_avg, volume)
        end if
 
        filename = trim(base_filename) // "_c.log"
 
        if (mype == 0 .and. .not. visited) then
        ! Delete the log when not doing a restart run
            if (restart_from_file == undefined .or. reset_time) then
                open(unit=my_unit,file=trim(filename),form='formatted',status='replace')
                if (mhd_bool == 1) then
                    write(my_unit,'(a)') 'global_time Tmax Tmin vmax B1max B2max B3max mag_avg'
                else
                    write(my_unit,'(a)') 'global_time Tmax Tmin vmax'
                end if
            end if
            visited = .true.
        end if
 
        if (mype == 0) then
            write(filename,"(a)") filename
            inquire(file=filename,exist=alive)
            if(alive) then
                open(unit=my_unit,file=filename,form='formatted',status='old',access='append')
            else
                open(unit=my_unit,file=filename,form='formatted',status='new')
            endif
 
            ! if number of output doubles is increase, don't forget to change the fmt_string above
            if (mhd_bool == 1) then
                write(my_unit, fmt_string) global_time, tmax, tmin, vmax, b1max, b2max, b3max, magn_avg
            else
                write(my_unit, fmt_string) global_time, tmax, tmin, vmax
            end if
            close(my_unit)
        end if
    end subroutine analytics_log
 
    pure function magnetic(w_vec, w_size) result(magn_energy)
        use mod_global_parameters
        integer, intent(in)          :: w_size
        double precision, intent(in) :: w_vec(w_size)
        double precision             :: magn_energy
 
        magn_energy = 0.5d0 * sum(w_vec(mag_mhd(:))**2,dim=ndir+1)
    end function magnetic
 
    ! Calculate both min and max of temperature on grid in one go.
    subroutine get_minmax_temperature(Tmax, Tmin)
        use mod_global_parameters
 
        double precision, intent(out) :: Tmax, Tmin
 
        integer                       :: iigrid, igrid
        double precision              :: Tmax_mype, Tmax_recv, Tmin_mype, Tmin_recv
        double precision              :: wlocal(ixG^T,1:nw), xlocal(ixG^T,1:ndim)
        double precision              :: Te(ixG^T), pth(ixG^T), rho(ixG^T), Rfactor(ixG^T)
 
        tmax_mype = -bigdouble
        tmin_mype = bigdouble
 
        !Loop over all the grids
        do iigrid = 1, igridstail
            igrid = igrids(iigrid)
 
            wlocal(ixg^t,1:nw) = ps(igrid)%w(ixg^t,1:nw)
            xlocal(ixg^t,1:ndim) = ps(igrid)%x(ixg^t,1:ndim)
            call eos%get_thermal_pressure(wlocal,xlocal,ixg^ll,ixg^ll,pth)
            call eos%get_rho(wlocal,xlocal,ixg^ll,ixm^ll,rho)
            call eos%get_Rfactor(wlocal,xlocal,ixg^ll,ixg^ll,rfactor)
            te(ixm^t) = pth(ixm^t)/(rho(ixm^t)*rfactor(ixm^t))
 
            ! Compare values on current grid to temporary max/min
            tmax_mype = max(tmax_mype,maxval(te(ixm^t)))
            tmin_mype = min(tmin_mype,minval(te(ixm^t)))
        end do
 
        ! Make the information available on all tasks
        call mpi_allreduce(tmax_mype, tmax_recv, 1, mpi_double_precision, &
            mpi_max, icomm, ierrmpi)
        call mpi_allreduce(tmin_mype, tmin_recv, 1, mpi_double_precision, &
            mpi_min, icomm, ierrmpi)
 
        tmax = tmax_recv
        tmin = tmin_recv
 
    end subroutine get_minmax_temperature
 
    subroutine get_max_velocity(vmax)
        use mod_global_parameters
        use mod_physics, only: phys_get_v
 
        double precision, intent(out) :: vmax
 
        integer                       :: iigrid, igrid
        double precision              :: vmax_mype,vmax_recv
        double precision              :: v_vec(ixG^T,1:ndir), v(ixG^T)
 
        vmax_mype = -bigdouble
 
        !Loop over all the grids
        do iigrid = 1, igridstail
            igrid = igrids(iigrid)
 
            call phys_get_v(ps(igrid)%w, ps(igrid)%x, ixg^ll, ixg^ll, v_vec)
 
            v(ixm^t) = sqrt(sum(v_vec(ixm^t,:)**2,dim=ndir+1))
 
            vmax_mype =max(vmax_mype ,maxval(v(ixm^t)))
 
        end do
 
        ! Make the information available on all tasks
        call mpi_allreduce(vmax_mype, vmax_recv, 1, mpi_double_precision, &
            mpi_max, icomm, ierrmpi)
 
        vmax = vmax_recv
 
    end subroutine get_max_velocity
 
    subroutine get_max_b(B1max, B2max, B3max)
        use mod_global_parameters
 
        double precision, intent(out) :: B1max, B2max, B3max
 
        integer                       :: iigrid, igrid
        double precision              :: B1max_mype, B1max_recv
        double precision              :: B2max_mype, B2max_recv
        double precision              :: B3max_mype, B3max_recv
 
        b1max_mype = -bigdouble
        b2max_mype = -bigdouble
        b3max_mype = -bigdouble
 
        !Loop over all the grids
        do iigrid = 1, igridstail
            igrid = igrids(iigrid)
 
            b1max_mype = max(b1max_mype, maxval(ps(igrid)%w(ixg^t, mag_mhd(1))))
            b2max_mype = max(b2max_mype, maxval(ps(igrid)%w(ixg^t, mag_mhd(2))))
            if (ndir >= 3) then
                b3max_mype = max(b3max_mype, maxval(ps(igrid)%w(ixg^t, mag_mhd(3))))
            end if
        end do
 
        ! Make the information available on all tasks
        call mpi_allreduce(b1max_mype, b1max_recv, 1, mpi_double_precision, &
            mpi_max, icomm, ierrmpi)
        call mpi_allreduce(b2max_mype, b2max_recv, 1, mpi_double_precision, &
            mpi_max, icomm, ierrmpi)
        call mpi_allreduce(b3max_mype, b3max_recv, 1, mpi_double_precision, &
            mpi_max, icomm, ierrmpi)
 
        b1max = b1max_recv
        b2max = b2max_recv
        b3max = b3max_recv
 
        if (b2max == -bigdouble) b2max = 0.d0
        if (b3max == -bigdouble) b3max = 0.d0
 
    end subroutine get_max_b
 
    subroutine symm(ixI^L,ixO^L,w,x,w_index,boundary)
    !> Symmetric boundary condition
        integer, intent(in) :: ixI^L, ixO^L, w_index
        double precision, intent(in) :: x(ixI^S, 1:ndim)
        double precision, intent(inout) :: w(ixI^S, 1:nw)
        character(len=*), intent(in) :: boundary
 
        select case (boundary)
        case ('outer')
            w(ixo^s, w_index) = w(ixomin1-1:ixomin1-nghostcells:-1^%1ixo^s, w_index)
        case ('inner')
            w(ixo^s, w_index) = w(ixomax1+nghostcells:ixomax1+1:-1^%1ixo^s, w_index)
        case default
            call mpistop('Unknown boundary side: ' // trim(boundary))
        end select
 
    end subroutine symm
 
    subroutine split_symm(ixI^L,ixO^L,w,x,w_index,equi_w,boundary)
    !> Symmetric boundary condition applied to the perturbation only
        integer, intent(in) :: ixI^L, ixO^L, w_index
        double precision, intent(in) :: x(ixI^S, 1:ndim)
        double precision, intent(inout) :: w(ixI^S, 1:nw)
        double precision, intent(in) :: equi_w(ixI^S)
        character(len=*), intent(in) :: boundary
 
        select case (boundary)
        case ('outer')
            w(ixo^s, w_index) = equi_w(ixo^s) &
                + (w(ixomin1-1:ixomin1-nghostcells:-1^%1ixo^s, w_index) - equi_w(ixomin1-1:ixomin1-nghostcells:-1^%1ixo^s))
        case ('inner')
            w(ixo^s, w_index) = equi_w(ixo^s) &
                + (w(ixomax1+nghostcells:ixomax1+1:-1^%1ixo^s, w_index) - equi_w(ixomax1+nghostcells:ixomax1+1:-1^%1ixo^s))
        case default
            call mpistop('Unknown boundary side: ' // trim(boundary))
        end select
 
    end subroutine split_symm
 
    subroutine extrapolate(ixI^L,ixO^L,w,x,w_index,boundary)
    !> Extrapolate boundary condition
        integer, intent(in) :: ixI^L, ixO^L, w_index
        double precision, intent(in) :: x(ixI^S, 1:ndim)
        double precision, intent(inout) :: w(ixI^S, 1:nw)
        character(len=*), intent(in) :: boundary
        integer :: ix1
 
        select case (boundary)
        case ('outer')
            do ix1 = ixomin1, ixomax1
                w(ix1^%1ixo^s, w_index) = 2.0d0 * w(ix1-1^%1ixo^s, w_index) &
                    - w(ix1-2^%1ixo^s, w_index)
            end do
        case ('inner')
            do ix1 = ixomax1, ixomin1, -1
                w(ix1^%1ixo^s, w_index) = 2.0d0 * w(ix1+1^%1ixo^s, w_index) &
                    - w(ix1+2^%1ixo^s, w_index)
            end do
        case default
            call mpistop('Unknown boundary side: ' // trim(boundary))
        end select
 
    end subroutine extrapolate
 
    subroutine split_extrapolate(ixI^L,ixO^L,w,x,w_index,equi_w,boundary)
    !> Extrapolate boundary condition applied to the perturbation only
        integer, intent(in) :: ixI^L, ixO^L, w_index
        double precision, intent(in) :: x(ixI^S, 1:ndim)
        double precision, intent(inout) :: w(ixI^S, 1:nw)
        double precision, intent(in) :: equi_w(ixI^S)
        character(len=*), intent(in) :: boundary
        integer :: ix1
 
        select case (boundary)
        case ('outer')
            do ix1 = ixomin1, ixomax1
                w(ix1^%1ixo^s, w_index) = equi_w(ix1^%1ixo^s) &
                    + 2.0d0 * (w(ix1-1^%1ixo^s, w_index) - equi_w(ix1-1^%1ixo^s)) &
                    - (w(ix1-2^%1ixo^s, w_index) - equi_w(ix1-2^%1ixo^s))
            end do
        case ('inner')
            do ix1 = ixomax1, ixomin1, -1
                w(ix1^%1ixo^s, w_index) = equi_w(ix1^%1ixo^s) &
                    + 2.0d0 * (w(ix1+1^%1ixo^s, w_index) - equi_w(ix1+1^%1ixo^s)) &
                    - (w(ix1+2^%1ixo^s, w_index) - equi_w(ix1+2^%1ixo^s))
            end do
        case default
            call mpistop('Unknown boundary side: ' // trim(boundary))
        end select
 
    end subroutine split_extrapolate
 
end module mod_gimli
!