module heat_flux_m !! Calculation of the heat flux according to the chosen model use MPI use netcdf use error_handling_grillix_m, only : handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER, GRILLIX_ERR_SOLVER3D use screen_io_m, only : get_stdout use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use parallel_map_m, only : parallel_map_t use solver_aligned3d_m, only : solve_aligned3d_adjoint, params_solver_aligned3d_t use variable_m, only : variable_t use inplane_operators_m, only : inplane_operators_t use equilibrium_storage_m, only : equilibrium_storage_t use constants_m, only : PI use boundaries_perp_m, only : set_perpbnds, extrapolate_ghost_points use boundaries_braginskii_m, only : boundaries_braginskii_t use penalisation_m, only : penalisation_t ! Parameters use facade_params_m, only : & tinfo_size use params_feature_selection_m, only: & mms_on use params_brag_model_m, only : & mass_ratio_ei, tratio use params_brag_pardiss_model_m, only : & chipar0_e, chipar0_i, heatflux_cutoff_e, heatflux_cutoff_i, & landau_numlorentzians_e, landau_numlorentzians_i, landau_flutter_lhs_on use params_brag_boundaries_parpen_m, only : & pen_epsinv, sheath_heattransfac_e, sheath_heattransfac_i use params_brag_parsolver_te_m, only : & parsolver_te_select, params_parsolver_te use params_brag_parsolver_ti_m, only : & parsolver_ti_select, params_parsolver_ti ! MMS solution and source for the heat flux use mms_braginskii_m, only : mms_sol_braginskii_qe, mms_source_braginskii_landau_e, & mms_sol_braginskii_qi, mms_source_braginskii_landau_i implicit none ! Defining the fit parameters for the landau heat flux model real(GP), dimension(1), private, target :: alpha_debug = & [0.01315_GP] real(GP), dimension(1), private, target :: beta_debug = & [0.2044_GP] ! alpha_debug and beta_debug are for testing and debugging purpose only, ! they do not represent a valid physical model real(GP), dimension(3), private, target :: alpha_3 = & [0.01315_GP, 0.924_GP, 14.1365_GP] real(GP), dimension(7), private, target :: alpha_7 = & [0.007438_GP, 0.6161_GP, 5.9804_GP, & 37.9822_GP, 234.3654_GP, 1466.4331_GP, 12981.4634_GP] real(GP), dimension(12), private, target :: alpha_12 = & [0.001424_GP, 0.20736_GP, 2.5653_GP, 14.927_GP, 79.305_GP, & 419.2399_GP, 2215.7233_GP, 11709.7857_GP, 61885.2763_GP, & 327392.6096_GP, 1773350.1566_GP, 16903628.3745_GP] real(GP), dimension(3), private, target :: beta_3 = & [0.2044_GP, 1.3587_GP, 8.9643_GP] real(GP), dimension(7), private, target :: beta_7 = & [0.1678_GP, 1.1106_GP, 5.6457_GP, & 33.1536_GP, 202.738_GP, 1254.2144_GP, 9275.3323_GP] real(GP), dimension(12), private, target :: beta_12 = & [0.09419_GP, 0.6741_GP, 2.9628_GP, 14.43958_GP, 75.1106_GP, & 395.8293_GP, 2090.8877_GP, 11049.1471_GP, 58392.0969_GP, & 308695.7371_GP, 1645460.1472_GP, 10794779.4293_GP] ! Two arrays for the initial guess for qcondt_temp (2D: lorentzians, points) real(GP), dimension(:,:), private, save, allocatable :: guess_e !! Initial guess for the electron solve real(GP), dimension(:,:), private, save, allocatable :: guess_i !! Initial guess for the ion solve public :: heat_flux_landau public :: heat_flux_braginskii public :: heat_flux_free_streaming contains subroutine heat_flux_landau(comm_handler, & equi, equi_on_cano, equi_on_stag, & mesh_cano, mesh_stag, & map, & penalisation_stag, & opsinplane_cano, opsinplane_stag, & boundaries, & stag_ne, stag_te, stag_ti, heat_gradpar_logte_stag, & heat_gradpar_logti_stag, tau, qcondt, kind_of_species, & apar, full_apar, & bnd_descrs_nmn_cano, iters, ress, true_ress) !! Contains a routine that calculates the heat flux !! according to the landau heat flux model !! Ref.: J.G. Chen, X.Q. Xu, Y.A. Lei, Extension of Landau-fluid closure !! to weakly collisional plasma regime, Computer Physics Communications, !! Volume 236, 2019, Pages 128-134, ISSN 0010-4655, !! https://doi.org/10.1016/j.cpc.2018.10.024. type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(inout) :: equi !! Equilibrium class(equilibrium_storage_t), intent(in) :: equi_on_cano !! Equilibrim quantities on canonical mesh class(equilibrium_storage_t), intent(in) :: equi_on_stag !! Equilibrim quantities on staggered mesh type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation_stag !! Penalisation (staggered) type(inplane_operators_t), intent(in) :: opsinplane_cano !! In-plane operators (canonical) type(inplane_operators_t), intent(in) :: opsinplane_stag !! In-plane operators (staggered) type(boundaries_braginskii_t), intent(in) :: boundaries !! Boundary information for the BRAGINSKII model real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: stag_ne !! Density on the staggered grid real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: stag_te !! Electron temperature on the staggered grid real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: stag_ti !! Ion temperature on the staggered grid real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: heat_gradpar_logte_stag !! Gradient of logarithm of electron temperature on the staggered grid real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: heat_gradpar_logti_stag !! Gradient of logarithm of ion temperature on the staggered grid real(GP), intent(in) :: tau !! Time real(GP), dimension(mesh_stag%get_n_points()), intent(inout) :: qcondt !! Heat flux character(len=4), intent(in) :: kind_of_species !! Kind of species for routine call ('elec' or 'ions') type(variable_t), intent(in) :: apar !! Parallel electromagnetic potential type(variable_t), intent(in) :: full_apar !! Parallel electromagnetic potential on the full grid integer, dimension(mesh_cano%get_n_points_boundary()), intent(in) :: bnd_descrs_nmn_cano !! Boundary descriptors for the helholtz solver with electromagnetic flutter integer, intent(out), dimension(12) :: iters !! Number of iterations of the Helmholtz solver (integers) real(GP), intent(out), dimension(12) :: ress !! Pseudo-residuals of the Helmholtz solver real(GP), intent(out), dimension(12) :: true_ress !! True residuals of the Helmholtz solver ! Local variables type(variable_t) :: stag_temperature ! Temperature of "kind_of_species" on the staggered grid type(variable_t) :: gradpar_log_temperature ! Logarithm of parallel gradient of the temperature for "kind_of_species" type(variable_t) :: qcondt_temp ! Temporary electron heat flux for the single lorentzians real(GP) :: qcondt_boundary ! Boundary condition for the heat flux in the helmholtz solver real(GP), dimension(:), pointer :: alpha ! Array with hard coded parameters alpha real(GP), dimension(:), pointer :: beta ! Array with hard coded parameters beta real(GP) :: delta ! Paramter to match the Braginskii limit for electrons real(GP) :: factor_l_zeta ! Factor that gets evaluated separatly to shorten equations real(GP) :: sqrt8bypi = sqrt(8.0_GP / PI) ! Factor sqrt(8/pi) integer :: number_of_lorentzians ! Number of lorentzians used in this call of the subroutine integer :: ki, kb, kg, l, n ! Running indices real(GP) :: pen_fac ! Penalisation factor real(GP) :: gamma_sheath ! Sheath factor integer, save :: cnt = 0 ! Counter for calls of the subroutine real(GP) :: factor_fix_norm ! Factor to correct the normalsation of this module ! Parameters for the parallel helmholtz solver real(GP), dimension(mesh_stag%get_n_points()) :: lambda_par real(GP), dimension(mesh_stag%get_n_points()) :: xi_par real(GP), dimension(mesh_stag%get_n_points()) :: rhs_par integer, dimension(mesh_stag%get_n_points_boundary()) :: bnd_descrs_qcondt integer :: cnt_matvec_parslv integer :: niter_parslv real(GP) :: res, true_res, flutter_fac_cano, flutter_fac_stag type(params_solver_aligned3d_t) :: params_parsolver character(len=:), allocatable :: parsolver_type ! Coordinates for MMS solution and source real(GP) :: x, y, phi_stag ! Function call counter cnt = cnt + 1 phi_stag = map%get_phi_stag() ! Allocate initial guess if (cnt == 1) then allocate(guess_e(landau_numlorentzians_e, mesh_stag%get_n_points())) allocate(guess_i(landau_numlorentzians_i, mesh_stag%get_n_points())) endif ! Initialize stag_temperature, gradpar_te and qcondte_temp and set them to the right grid call stag_temperature%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.) call gradpar_log_temperature%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.) call qcondt_temp%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.) ! Define all variables that are different for the electron and ion cases if (kind_of_species == 'elec') then delta = 0.50_GP factor_fix_norm = 1.0_GP / sqrt(mass_ratio_ei) factor_l_zeta = 3.16_GP * delta / (chipar0_e * sqrt(mass_ratio_ei)) number_of_lorentzians = landau_numlorentzians_e params_parsolver = params_parsolver_te gamma_sheath = sheath_heattransfac_e ! Copy solver parameters parsolver_type = parsolver_te_select params_parsolver%krylov_method = params_parsolver_te%krylov_method params_parsolver%precond_type = params_parsolver_te%precond_type params_parsolver%resmax = params_parsolver_te%resmax params_parsolver%maxiter = params_parsolver_te%maxiter params_parsolver%nrestart = params_parsolver_te%nrestart params_parsolver%dbgout = params_parsolver_te%dbgout else if (kind_of_species == 'ions') then delta = 0.41_GP factor_fix_norm = sqrt(tratio) factor_l_zeta = 3.9_GP * delta * sqrt(tratio) / chipar0_i number_of_lorentzians = landau_numlorentzians_i params_parsolver = params_parsolver_ti gamma_sheath = sheath_heattransfac_i ! Copy solver parameters parsolver_type = parsolver_ti_select params_parsolver%krylov_method = params_parsolver_ti%krylov_method params_parsolver%precond_type = params_parsolver_ti%precond_type params_parsolver%resmax = params_parsolver_ti%resmax params_parsolver%maxiter = params_parsolver_ti%maxiter params_parsolver%nrestart = params_parsolver_ti%nrestart params_parsolver%dbgout = params_parsolver_ti%dbgout else call handle_error('Only ions or elec valid as species', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t(kind_of_species)) endif ! Initialize the vectors of values for alpha and beta ! (fit parameters for the lorentzian functions) if (number_of_lorentzians == 1) then alpha => alpha_debug beta => beta_debug elseif (number_of_lorentzians == 3) then alpha => alpha_3 beta => beta_3 else if (number_of_lorentzians == 7) then alpha => alpha_7 beta => beta_7 else if (number_of_lorentzians == 12) then alpha => alpha_12 beta => beta_12 else call handle_error('Number of Lorentzians invalid, only (3,7,12) allowed', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Nlorentzians ', & [number_of_lorentzians])) endif !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, kind_of_species, & !$omp stag_te, stag_ti, stag_temperature, gradpar_log_temperature, & !$omp heat_gradpar_logte_stag, heat_gradpar_logti_stag, qcondt, qcondt_temp) !$omp do do l = 1, mesh_stag%get_n_points() ! Define field quantities that are different for the electron and ion cases if (kind_of_species == 'elec') then stag_temperature%vals(l) = stag_te(l) gradpar_log_temperature%vals(l) = heat_gradpar_logte_stag(l) else if (kind_of_species == 'ions') then stag_temperature%vals(l) = stag_ti(l) gradpar_log_temperature%vals(l) = heat_gradpar_logti_stag(l) endif ! Initialize the heat conduction array, because we add upon this ! and the temporary one, so that it is well defined for the first solve qcondt(l) = 0.0_GP qcondt_temp%vals(l) = 0.0_GP enddo !$omp enddo !$omp end parallel ! Loop over the single lorentzians do n = 1, number_of_lorentzians ! Calculate solver parameters for interior !$omp parallel default(none) & !$omp private(ki, kb, kg, l, qcondt_boundary, pen_fac) & !$omp shared(cnt, n, mesh_stag, penalisation_stag, bnd_descrs_qcondt, & !$omp beta, tratio, gamma_sheath, pen_epsinv, & !$omp kind_of_species, alpha, sqrt8bypi, factor_l_zeta, factor_fix_norm, & !$omp lambda_par, xi_par, rhs_par, guess_e, guess_i, qcondt_temp, & !$omp stag_ne, stag_te, stag_ti, stag_temperature, gradpar_log_temperature, & !$omp boundaries) !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) ! Select the initial guess for qcondt_temp (zero at first solve) if (cnt > 2) then if (kind_of_species == 'elec') then qcondt_temp%vals(l) = guess_e(n,l) elseif (kind_of_species == 'ions') then qcondt_temp%vals(l) = guess_i(n,l) endif endif if (n == 1) then qcondt_boundary = gamma_sheath * stag_ne(l) * stag_temperature%vals(l) * & sqrt(stag_te(l) + tratio * stag_ti(l)) * & penalisation_stag%get_dirindfun_val(ki) / & (factor_l_zeta * alpha(n) * sqrt8bypi) else qcondt_boundary = 0.0_GP endif pen_fac = 1.0_GP - penalisation_stag%get_charfun_val(ki) lambda_par(l) = (pen_fac * (factor_l_zeta * (stag_ne(l) / & stag_temperature%vals(l)**2) * beta(n))**2 + & pen_epsinv * & penalisation_stag%get_charfun_val(ki)) / & (factor_l_zeta * alpha(n) * sqrt8bypi) xi_par(l) = pen_fac / (factor_l_zeta * alpha(n) * sqrt8bypi) rhs_par(l) = - pen_fac * factor_fix_norm * stag_ne(l)**2 / & sqrt(stag_temperature%vals(l)) * & gradpar_log_temperature%vals(l) + & pen_epsinv * & penalisation_stag%get_charfun_val(ki) * & qcondt_boundary enddo !$omp end do ! Calculate solver parameters for boundary !$omp do do kb = 1, mesh_stag%get_n_points_boundary() l = mesh_stag%boundary_indices(kb) lambda_par(l) = 1.0_GP xi_par(l) = 0.0_GP if (kind_of_species == 'elec') then bnd_descrs_qcondt(kb) = boundaries%qe%types(kb) rhs_par(l) = boundaries%qe%vals(kb) ! Select electron initial guess for qcondt_temp (zero at first solve) if (cnt > 2) then qcondt_temp%vals(l) = guess_e(n,l) endif elseif (kind_of_species == 'ions') then bnd_descrs_qcondt(kb) = boundaries%qi%types(kb) rhs_par(l) = boundaries%qi%vals(kb) ! Select ion initial guess for qcondt_temp (zero at first solve) if (cnt > 2) then qcondt_temp%vals(l) = guess_i(n,l) endif endif enddo !$omp end do ! Calculate solver parameters for ghost cells !$omp do do kg = 1, mesh_stag%get_n_points_ghost() l = mesh_stag%ghost_indices(kg) ! Select the initial guess for qcondt_temp (zero at first solve) if (cnt > 2) then if (kind_of_species == 'elec') then qcondt_temp%vals(l) = guess_e(n,l) elseif (kind_of_species == 'ions') then qcondt_temp%vals(l) = guess_i(n,l) endif endif lambda_par(l) = 1.0_GP xi_par(l) = 0.0_GP rhs_par(l) = 0.0_GP enddo !$omp end do !$omp end parallel if (mms_on) then if (kind_of_species == 'elec') then !$omp parallel default(none) & !$omp private(ki, kb, kg, l, x, y) & !$omp shared(equi, tau, mesh_stag, phi_stag, & !$omp n, alpha, beta, bnd_descrs_qcondt, rhs_par, & !$omp boundaries) !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) rhs_par(l) = rhs_par(l) + mms_source_braginskii_landau_e(equi, & x, y, phi_stag, tau, alpha(n), beta(n)) enddo !$omp end do !$omp do do kb = 1, mesh_stag%get_n_points_boundary() l = mesh_stag%boundary_indices(kb) bnd_descrs_qcondt(kb) = boundaries%qe%types(kb) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) rhs_par(l) = mms_sol_braginskii_qe(equi, x, y, phi_stag, tau) enddo !$omp end do !$omp do do kg = 1, mesh_stag%get_n_points_ghost() l = mesh_stag%ghost_indices(kg) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) rhs_par(l) = mms_sol_braginskii_qe(equi, x, y, phi_stag, tau) enddo !$omp end do !$omp end parallel else if (kind_of_species == 'ions') then !$omp parallel private(ki, kb, kg, l, x, y) & !$omp shared(equi, tau, mesh_stag, phi_stag, & !$omp n, alpha, beta, bnd_descrs_qcondt, rhs_par, & !$omp boundaries) !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) rhs_par(l) = rhs_par(l) + mms_source_braginskii_landau_i(equi, & x, y, phi_stag, tau, alpha(n), beta(n)) enddo !$omp end do !$omp do do kb = 1, mesh_stag%get_n_points_boundary() l = mesh_stag%boundary_indices(kb) bnd_descrs_qcondt(kb) = boundaries%qi%types(kb) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) rhs_par(l) = mms_sol_braginskii_qi(equi, x, y, phi_stag, tau) enddo !$omp end do !$omp do do kg = 1, mesh_stag%get_n_points_ghost() l = mesh_stag%ghost_indices(kg) x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) rhs_par(l) = mms_sol_braginskii_qi(equi, x, y, phi_stag, tau) enddo !$omp end do !$omp end parallel endif endif ! Solve the helmholtz problem if (landau_flutter_lhs_on) then flutter_fac_cano = 1.0_GP flutter_fac_stag = 1.0_GP else flutter_fac_cano = 0.0_GP flutter_fac_stag = 0.0_GP endif call solve_aligned3d_adjoint(comm_handler, & equi_on_cano=equi_on_cano, equi_on_stag=equi_on_stag, & mesh_cano=mesh_cano, mesh_stag=mesh_stag, & map=map, & params_solver=params_parsolver, solver_type=parsolver_type, & rhs=rhs_par, sol=qcondt_temp%vals, & niter=niter_parslv, res=res, true_res=true_res,cnt_matvec=cnt_matvec_parslv, & lambda=lambda_par, xi=xi_par, & bnd_descrs=bnd_descrs_qcondt, & flutter_fac_cano=flutter_fac_cano, flutter_fac_stag=flutter_fac_stag, & opsinplane_cano=opsinplane_cano, opsinplane_stag=opsinplane_stag, & apar_stag=apar%vals, apar_cano=full_apar%vals, & bnd_descrs_flux=bnd_descrs_nmn_cano) !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, kind_of_species, n, & !$omp qcondt, qcondt_temp, guess_e, guess_i) !$omp do do l = 1, mesh_stag%get_n_points() ! Sum up the individual contributions of the ! Lorentzians for the total heat flux qcondt(l) = qcondt(l) + qcondt_temp%vals(l) ! Store the value of qcondt_temp as initial guess for the next timestep if (kind_of_species == 'elec') then guess_e(n,l) = qcondt_temp%vals(l) elseif (kind_of_species == 'ions') then guess_i(n,l) = qcondt_temp%vals(l) endif enddo !$omp end do !$omp end parallel ! Write the solver values to tinfo if (niter_parslv > 0) then iters(n) = cnt_matvec_parslv else iters(n) = niter_parslv endif ! Write residuum to rinfo ress(n) = res true_ress(n) = true_res enddo if (mms_on) then ! The solution for the heatflux needs to be divided by the number of lorentzians ! since each lorentzian tries to find the MMS solution for qcondt !$omp parallel private(l) !$omp do do l = 1, mesh_stag%get_n_points() qcondt(l) = qcondt(l) / number_of_lorentzians enddo !$omp end do !$omp end parallel endif !$omp parallel default(none) shared(mesh_stag, qcondt) call extrapolate_ghost_points(mesh_stag, qcondt) !$omp end parallel end subroutine subroutine heat_flux_braginskii(mesh_stag, stag_temperature, fac_free_streaming, & gradpar_log_temperature, stag_ne, qcondt, kind_of_species, bnd_descrs_stag) !! Computes heat flux according to Braginskii model with free streaming limter type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: stag_temperature !! Temperatures on the staggered grid real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: stag_ne !! Density on the staggered grid real(GP), intent(in) :: fac_free_streaming !! Calculated free streaming factor real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: gradpar_log_temperature !! Gradients of the logarithms of temperature real(GP), dimension(mesh_stag%get_n_points()), intent(inout) :: qcondt !! Heat fluxes as Output character(len=4), intent(in) :: kind_of_species !! Kind of species for routine call ('elec' or 'ions') integer, dimension(mesh_stag%get_n_points_boundary()), intent(in) :: bnd_descrs_stag !! Boundary condition for heat flux real(GP) :: heatflux_coeff_e, heatflux_coeff_i integer :: l if (kind_of_species == 'elec') then !$omp parallel default(none) private(l, heatflux_coeff_e) & !$omp shared(mesh_stag, bnd_descrs_stag, & !$omp chipar0_e, fac_free_streaming, heatflux_cutoff_e, & !$omp stag_ne, stag_temperature, gradpar_log_temperature, qcondt) !$omp do do l = 1, mesh_stag%get_n_points() ! Parallel heat flux coefficient (electrons) ! Braginskii heat flux heatflux_coeff_e = chipar0_e * sqrt(stag_temperature(l)**5) ! Free streaming limiter heatflux_coeff_e = heatflux_coeff_e / (1.0_GP + fac_free_streaming * & stag_temperature(l)**2 / stag_ne(l)) ! Hard heat flux limiter heatflux_coeff_e = min(heatflux_coeff_e, heatflux_cutoff_e) qcondt(l) = - heatflux_coeff_e * stag_temperature(l) * & gradpar_log_temperature(l) enddo !$omp end do call set_perpbnds(mesh_stag, bnd_descrs_stag, qcondt) call extrapolate_ghost_points(mesh_stag, qcondt) !$omp end parallel elseif (kind_of_species == 'ions') then !$omp parallel default(none) private(l, heatflux_coeff_i) & !$omp shared(mesh_stag, bnd_descrs_stag, & !$omp chipar0_i, heatflux_cutoff_i, fac_free_streaming, & !$omp stag_ne, stag_temperature, gradpar_log_temperature, qcondt) !$omp do do l = 1, mesh_stag%get_n_points() ! Parallel heat flux coefficient (ions) ! Braginskii heat flux heatflux_coeff_i = chipar0_i * sqrt(stag_temperature(l)**5) ! Free streaming limiter heatflux_coeff_i = heatflux_coeff_i / (1.0_GP + fac_free_streaming * & stag_temperature(l)**2 / stag_ne(l)) ! Hard heat flux limiter heatflux_coeff_i = min(heatflux_coeff_i, heatflux_cutoff_i) qcondt(l) = - heatflux_coeff_i * stag_temperature(l) * & gradpar_log_temperature(l) enddo !$omp end do call set_perpbnds(mesh_stag, bnd_descrs_stag, qcondt) call extrapolate_ghost_points(mesh_stag, qcondt) !$omp end parallel endif end subroutine subroutine heat_flux_free_streaming(mesh_stag, stag_temperature, stag_ne, & gradpar_log_temperature, qcondt, kind_of_species, bnd_descrs_stag) !! Computes heat flux according to free streaming model type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: stag_temperature !! Temperatures on the staggered grid real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: stag_ne !! Density on the staggered grid real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: gradpar_log_temperature !! Gradients of the logarithms of temperature real(GP), dimension(mesh_stag%get_n_points()), intent(inout) :: qcondt !! Heat fluxes as Output character(len=4), intent(in) :: kind_of_species !! Kind of species for routine call ('elec' or 'ions') integer, dimension(mesh_stag%get_n_points_boundary()), intent(in) :: bnd_descrs_stag !! Boundary condition for heat flux real(GP) :: heatflux_coeff_e, heatflux_coeff_i integer :: l if (kind_of_species == 'elec') then !$omp parallel default(none) private(l, heatflux_coeff_e) & !$omp shared(mesh_stag, bnd_descrs_stag, mass_ratio_ei, & !$omp stag_ne, stag_temperature, gradpar_log_temperature, qcondt) !$omp do do l = 1, mesh_stag%get_n_points() ! Parallel heat flux coefficient (electrons) heatflux_coeff_e = stag_ne(l) * sqrt(stag_temperature(l)) / sqrt(mass_ratio_ei) qcondt(l) = - heatflux_coeff_e * stag_temperature(l) * gradpar_log_temperature(l) enddo !$omp end do call set_perpbnds(mesh_stag, bnd_descrs_stag, qcondt) call extrapolate_ghost_points(mesh_stag, qcondt) !$omp end parallel elseif(kind_of_species == 'ions') then !$omp parallel private(l, heatflux_coeff_i) & !$omp shared(mesh_stag, bnd_descrs_stag, tratio, & !$omp stag_ne, stag_temperature, gradpar_log_temperature, qcondt) !$omp do do l = 1, mesh_stag%get_n_points() ! Parallel heat flux coefficient (ions) heatflux_coeff_i = stag_ne(l) * sqrt(stag_temperature(l)) * sqrt(tratio) qcondt(l) = - heatflux_coeff_i * stag_temperature(l) * gradpar_log_temperature(l) enddo !$omp end do call set_perpbnds(mesh_stag, bnd_descrs_stag, qcondt) call extrapolate_ghost_points(mesh_stag, qcondt) !$omp end parallel endif end subroutine end module