hyperdiffusion.f90

Go to the documentation of this file.

module hyperdiffusion

  use ids_schemas, only: ids_real, ids_real_invalid
  use type_transport_solver_numerics, only: hyperdiffusion_type

contains

  subroutine get_hyper_values(hyper, rho, nrho_bnd, prof, deriv, diff, params, active, prof_name, error_index, error_message)
    use ets_plasma, only: type_hyperdiffusion_params
    use control_parameters, only: id_diffusion_and_advection, id_diffusion_and_source

    ! Input/Output
    type(hyperdiffusion_type),        intent(inout) :: hyper
    real(ids_real),                   intent(in)    :: rho(:)
    integer,                          intent(in)    :: nrho_bnd
    real(ids_real),                   intent(in)    :: prof(:)
    real(ids_real),                   intent(in)    :: deriv(:)
    real(ids_real),                   intent(in)    :: diff(:)
    type(type_hyperdiffusion_params), intent(in)    :: params
    logical,                          intent(in)    :: active
    character(*),                     intent(in)    :: prof_name
    integer,                          intent(out)   :: error_index
    character(*),                     intent(inout) :: error_message

    ! Internal
    real(ids_real) :: max_rho(size(rho))
    integer :: nrho

    nrho = size(rho)

    hyper%diff  = 0.0_ids_real
    hyper%vconv = 0.0_ids_real
    hyper%src   = 0.0_ids_real

    ! If hyper-diffusion is turned off, then we should return zero coefficients
    if (.not. active) return

    ! Set the value of the hyper diffusion, but do not fill the entire rho-domain...
    ! - leave at least one grid point by boundary between the core and the edge, here described by nrho_bnd;
    ! - do not fill the hyper-diffusion beyond the "rho_norm_cutoff".
    max_rho = min( 0.5_ids_real*(rho(nrho_bnd-1)+rho(nrho_bnd)) , params%rho_norm_cutoff)
    where (rho <= max_rho)
        hyper%diff = params%explicit + params%implicit * maxval(diff(1:nrho_bnd))
    end where

    ! Fill the term that will cancel the flux from the hyper diffusion, i.e. either a convection term or an source term
    if (params%terms == id_diffusion_and_advection) THEN
        where (abs(deriv) <= tiny(1.0_ids_real))
            hyper%vconv = hyper%diff * deriv / prof
        end where
    elseif (params%terms == id_diffusion_and_source) THEN
        error_index   = -1
        error_message = 'in hyper_value ('//prof_name//'), hyperdiffusion with sources is not yet implemented'
        return
    else
        error_index   = -1
        error_message = 'in hyper_value ('//prof_name//'), cannot determine how to set up the hyperdiffusion'
        return
    end if

    error_index   = 0
 
  end subroutine get_hyper_values

end module hyperdiffusion