timestep_zhdanov_m.f90 Source File

module timestep_zhdanov_m
    !! Implementation of ZHDANOV model
    use comm_handler_m,             only : comm_handler_t
    use precision_grillix_m,        only : GP, GP_NAN
    use parallel_map_m,             only : parallel_map_t
    use penalisation_m,             only : penalisation_t
    use polars_m,                   only : polars_t
    use equilibrium_storage_m,      only : equilibrium_storage_t
    use variable_m,                 only : variable_t
    use multistep_m,                only : multistep_storage_t, karniadakis_t
    use compute_electron_density_m, only : compute_electron_density 
    use params_zhdanov_general_m,   only : params_zhdanov_general_t
    use params_zhdanov_species_m,   only : params_zhdanov_species_t
    use equilibrium_m,              only : equilibrium_t
    use mesh_cart_m,                only : mesh_cart_t  
    use params_zhdanov_general_m,   only : params_zhdanov_general_t
    use params_zhdanov_species_m,   only : params_zhdanov_species_t
    use mms_zhdanov_m,              only : MmsSourceDens
    implicit none

    public :: timestep_zhdanov

    ! Auxiliary variables
    type(variable_t), private, save :: log_dens_species
    !! Logarithmic density

    integer, private, save :: cnt = 0
    !! Counter for calling timestep_zhdanov

contains

    subroutine timestep_zhdanov(comm_handler, equi, equi_storage, mesh,&
                                   map, penalisation, polars,             &
                                   params_general,                        &
                                   params_species,                        &
                                   tau,                                   &
                                   dens_species,                          &
                                   tstep_continuity)
        !! Evolves variable a single time-step according to zhdanov model
        type(comm_handler_t), intent(in) :: comm_handler
        !! Communicators
        class(equilibrium_t), intent(inout) :: equi
        !! Equilibrium
        type(equilibrium_storage_t), intent(in) :: equi_storage
        !! Equilibrium storage enabling faster performance at certain locations
        type(mesh_cart_t), intent(inout) :: mesh
        !! Mesh within poloidal plane
        type(parallel_map_t), intent(in) :: map
        !! Parallel map
        type(penalisation_t), intent(in) :: penalisation
        !! Penalisation    
        type(polars_t), intent(in) :: polars
        !! Polar grid and operators
        type(params_zhdanov_general_t), intent(in)  :: params_general
        !! General parameters related with ZHDANOV model
        type(params_zhdanov_species_t), intent(in)  :: params_species
        !! Species related parameters related with ZHDANOV model
        real(GP), intent(in) :: tau
        !! Time
        type(variable_t), intent(inout):: dens_species
        !! Species density   
        type(karniadakis_t), intent(inout) :: tstep_continuity
        !! Time-step integrator for ion density

        real(GP), dimension(mesh%get_n_points()) :: d_log_dens_species
        !! Change of variable, i.e. rhs (explicit part)) of partial differential equation

        real(GP), dimension(mesh%get_n_points()) :: dens_species_adv, &
        log_dens_species_adv 
        !! Variables advanced in time at t+1
        real(GP), dimension(mesh%get_n_points()) :: charge_dens
        !! Z-weighted species density
                
        integer :: l
        real(GP) :: phi, x, y

        cnt = cnt + 1        

        phi = mesh%get_phi()

        ! Setup structures at first call ----------------------------------------------------------
        if (cnt == 1) then
            call setup_work_variables()

            ! Initialise logarithmic quantities (for subsequent time steps, they are pre-computed at the end of the routine).
            !$omp parallel do 
            do l = 1, mesh%get_n_points()
                log_dens_species%vals(l) = log( dens_species%vals(l) )
            enddo
            !$omp end parallel do
        endif
      
        if (comm_handler%get_species() /= 0)  then
            ! Set profiles dens_species
            !$omp parallel private(l)
            !$omp do
            do l = 1, mesh%get_n_points()
                ! Continuity equation (explicit part) -----------------------------
                d_log_dens_species(l) = params_species%source_param / &
                   dens_species%vals(l)
            enddo
            !$omp end do
            !$omp end parallel
            
            if (params_general%mms_on) then
                !$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)
                    d_log_dens_species(l) = d_log_dens_species(l) + &
                        MmsSourceDens(equi, comm_handler%get_species(), x, y,&
                        phi, tau, params_species) / dens_species%vals(l)
                enddo
                !$omp end do
                !$omp end parallel
            endif
           
            ! Explicit time advancement of variables ------------------------------
            call tstep_continuity%advance(d_log_dens_species, &
                                          log_dens_species%vals, &
                                          log_dens_species_adv)
           
            call tstep_continuity%shift_storage((/log_dens_species%vals, &
                                                d_log_dens_species/))

            !$omp parallel private(l)
            ! Compute actual species densities from their logarithms
            !$omp do 
            do l = 1, mesh%get_n_points()
                dens_species_adv(l) = exp( log_dens_species_adv(l) )
            enddo
            !$omp end do 
            !$omp end parallel

        else
            ! Fill electron density with zero 
            ! -> will be obtained via compute_electron_density
            !$omp parallel do private(l)
            do l = 1, mesh%get_n_points()
                dens_species_adv(l) = 0.0_GP
            enddo
            !$omp end parallel do
                                            
        endif

        call compute_electron_density(comm_handler, mesh, &
                                      params_species, dens_species_adv) 

        !$omp parallel do private(l)
        do l = 1, mesh%get_n_points()
            dens_species%vals(l) = dens_species_adv(l)
            ! These fields are saved within the timestep module for the next time step.
            log_dens_species%vals(l) = log_dens_species_adv(l)
        enddo
        !$omp end parallel do

    contains
        
        subroutine setup_work_variables()
            !! Allocates work variables
            call log_dens_species%init(comm_handler, mesh%get_n_points(), &
               staggered = .false.)
        end subroutine

    end subroutine

end module