module penalisation_neutrals_m ! Module for calculating penalisation values for the neutrals module use perf_m, only : perf_start, perf_stop use error_handling_grillix_m, only : handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use descriptors_neutrals_m, only : BND_TYPE_DIRICHLET_ZERO, & BND_TYPE_NEUMANN, & BND_TYPE_NONE use descriptors_m, only : DISTRICT_CORE, DISTRICT_WALL, DISTRICT_DOME, DISTRICT_OUT use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use parallel_map_m, only : parallel_map_t use penalisation_m, only : penalisation_t use variable_m, only: variable_t use precision_grillix_m, only : GP, GP_NAN, GP_EPS use screen_io_m, only : get_stdout ! Parameters use params_neut_boundaries_parpen_m, only : second_pen_flag, & bnddescr_neut_dens_pen, bndval_neut_dens_peninto, bndval_neut_dens_penout, & bndval_neut_dens_peninto_second, bndval_neut_dens_penout_second, & bnddescr_neut_parmom_pen, bndval_neut_parmom_peninto, bndval_neut_parmom_penout, & bndval_neut_parmom_peninto_second, bndval_neut_parmom_penout_second, & bnddescr_neut_pressure_pen, bndval_neut_pressure_peninto, bndval_neut_pressure_penout, & bndval_neut_pressure_peninto_second, bndval_neut_pressure_penout_second, & on_neutrals_temp implicit none private public :: neutrals_dens_penalisation public :: neutrals_parmom_penalisation public :: neutrals_pressure_penalisation private :: apply_dirichlet_neutrals private :: apply_neumann_neutrals private :: apply_none contains subroutine neutrals_dens_penalisation(equi, mesh_cano, map, penalisation_cano, & neutrals_dens, neutrals_dens_pen_vals) !! Sets the penalisation values for neutrals density class(equilibrium_t), intent(in) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation_cano !! Penalisation (canonical) type(variable_t), intent(in) :: neutrals_dens !! Neutrals density real(GP), intent(out), dimension(mesh_cano%get_n_points_inner()) :: neutrals_dens_pen_vals !! Penalisation values for the neutrals density select case(bnddescr_neut_dens_pen) case (BND_TYPE_DIRICHLET_ZERO) call apply_dirichlet_neutrals(mesh_cano, equi, penalisation_cano, bndval_neut_dens_peninto, & bndval_neut_dens_penout, bndval_neut_dens_peninto_second, & bndval_neut_dens_penout_second, neutrals_dens_pen_vals) case (BND_TYPE_NEUMANN) call apply_neumann_neutrals(.false., mesh_cano, map, penalisation_cano, & neutrals_dens, neutrals_dens_pen_vals) case (BND_TYPE_NONE) call apply_none(mesh_cano, neutrals_dens, neutrals_dens_pen_vals) case default call handle_error('Pen_bndtype not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Boundary_condition_type: ', & [bnddescr_neut_dens_pen])) end select end subroutine subroutine neutrals_parmom_penalisation(equi, mesh_stag, map, penalisation_stag, & neutrals_parmom, neutrals_parmom_pen_vals) !! Sets the penalisation values for neutrals parallel momentum class(equilibrium_t), intent(in) :: equi !! Equilibrium 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(variable_t), intent(in) :: neutrals_parmom !! Neutrals parallel momentum real(GP), intent(out), dimension(mesh_stag%get_n_points_inner()) :: neutrals_parmom_pen_vals !! Penalisation values for the neutrals parallel momentum select case(bnddescr_neut_parmom_pen) case (BND_TYPE_DIRICHLET_ZERO) call apply_dirichlet_neutrals(mesh_stag, equi, penalisation_stag, bndval_neut_parmom_peninto, & bndval_neut_parmom_penout, bndval_neut_parmom_peninto_second, & bndval_neut_parmom_penout_second, neutrals_parmom_pen_vals) case default call handle_error('Pen_bndtype not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Boundary_condition_type: ', & [bnddescr_neut_parmom_pen])) end select end subroutine subroutine neutrals_pressure_penalisation(equi, mesh_cano, map, penalisation_cano, neutrals_dens, & neutrals_temp, neutrals_pressure, neutrals_pressure_pen_vals) !! Sets the penalisation values for neutrals pressure class(equilibrium_t), intent(in) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation_cano !! Penalisation (canonical) type(variable_t), intent(in) :: neutrals_dens !! Neutrals density type(variable_t), intent(in) :: neutrals_temp !! Neutrals temperature type(variable_t), intent(in) :: neutrals_pressure !! Neutrals pressure real(GP), intent(out), dimension(mesh_cano%get_n_points_inner()) :: neutrals_pressure_pen_vals !! Penalisation values for the neutrals pressure integer :: l, ki select case(bnddescr_neut_pressure_pen) case (BND_TYPE_DIRICHLET_ZERO) call apply_dirichlet_neutrals(mesh_cano, equi, penalisation_cano, bndval_neut_pressure_peninto, & bndval_neut_pressure_penout, bndval_neut_pressure_peninto_second, & bndval_neut_pressure_penout_second, neutrals_pressure_pen_vals) case (BND_TYPE_NEUMANN) if (on_neutrals_temp) then ! Apply neumann condition on temperature instead call apply_neumann_neutrals(.false., mesh_cano, map, penalisation_cano, & neutrals_temp, neutrals_pressure_pen_vals) ! Convert temperature back to pressure !$omp parallel do default(none) private(l, ki) & !$omp shared(mesh_cano, neutrals_pressure_pen_vals, neutrals_dens) do ki = 1, mesh_cano%get_n_points_inner() l = mesh_cano%inner_indices(ki) neutrals_pressure_pen_vals(ki) = neutrals_pressure_pen_vals(ki) * neutrals_dens%vals(l) enddo !$omp end parallel do else call apply_neumann_neutrals(.false., mesh_cano, map, penalisation_cano, & neutrals_pressure, neutrals_pressure_pen_vals) endif case default call handle_error('Pen_bndtype not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Boundary_condition_type: ', & [bnddescr_neut_pressure_pen])) end select end subroutine subroutine apply_dirichlet_neutrals(mesh_any, equi, penalisation_any, & pen_bndval_into, pen_bndval_out, pen_bndval_into_second, & pen_bndval_out_second, penalised_values) !! Calculates values for a zeroth order dirichlet boundary condition (i.e. direct value set) !! Penalisation values obtained from linear combination such that the following rules hold: !! pen_dirind = 1.0 ==> penalised_value = pen_bndval_into !! pen_dirind = -1.0 ==> penalised_value = pen_bndval_out !! pen_dirind = 0.0 ==> penalised_value = (pen_bndval_into + pen_bndval_out)/2 type(mesh_cart_t), intent(in) :: mesh_any !! Mesh (cano or stag) class(equilibrium_t), intent(in) :: equi !! Equilibrium type(penalisation_t), intent(in) :: penalisation_any !! Penalisation (consistent with mesh) real(GP), intent(in) :: pen_bndval_into !! Value or scale factor for penalisation boundary where magnetic field is directed towards target real(GP), intent(in) :: pen_bndval_out !! Value or scale factor for penalisation boundary where magnetic field is directed out of target real(GP), intent(in) :: pen_bndval_into_second !! Value or scale factor for secondary penalisation boundary where magnetic field is directed towards target real(GP), intent(in) :: pen_bndval_out_second !! Value or scale factor for secondary penalisation boundary where magnetic field is directed out of target real(GP), intent(out), dimension(mesh_any%get_n_points_inner()) :: penalised_values !! Array with dirichlet-zero values integer :: ki real(GP) :: pen_dirind !$omp parallel default(none) private(ki, pen_dirind) & !$omp shared(mesh_any, equi, penalisation_any, pen_bndval_into, pen_bndval_out, & !$omp pen_bndval_into_second, pen_bndval_out_second, penalised_values, second_pen_flag) !$omp do do ki = 1, mesh_any%get_n_points_inner() pen_dirind = penalisation_any%get_dirindfun_val(ki) penalised_values(ki) = pen_bndval_into * (1.0_GP + pen_dirind) / 2.0_GP & + pen_bndval_out * (1.0_GP - pen_dirind) / 2.0_GP if (second_pen_flag /= 0) then ! Apply penalised values to the secondary penalisation regions if ((second_pen_flag * (mesh_any%get_y(ki) - equi%y0)) > 0.0_GP) then ! If second_pen_flag= 1, apply to the region with y > y0 ! If second_pen_flag=-1, apply to the region with y < y0 penalised_values(ki) = pen_bndval_into_second * (1.0_GP + pen_dirind) / 2.0_GP & + pen_bndval_out_second * (1.0_GP - pen_dirind) / 2.0_GP endif endif enddo !$omp end do !$omp end parallel end subroutine subroutine apply_neumann_neutrals(staggered, mesh, map, penalisation, values, penalised_values, gradients) !! Calculate penalisation values for a local neumann condition, by taking the upstream point value logical, intent(in) :: staggered !! Whether work variable is defined on staggered mesh, otherwise on canonical type(mesh_cart_t), intent(in) :: mesh !! Mesh (any, consistent with staggering) type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation !! Penalisation (any, consistent with staggering) type(variable_t), intent(in) :: values !! Variable with filled halos real(GP), optional, dimension(mesh%get_n_points_inner()) :: gradients !! Point-wise gradients, defined only on the grid (indexed by ka) real(GP), intent(out), dimension(mesh%get_n_points_inner()) :: penalised_values !! Array with neumann-zero values, defined only on the grid (indexed by ki) integer :: ki, l real(GP) :: pen_dirind real(GP), dimension(mesh%get_n_points()) :: vals_fwd, vals_bwd real(GP) :: grad, fwd_bwd_average real(GP) :: map_dpar_fwd, map_dpar_bwd !$omp parallel default(none) & !$omp private(ki, l, grad, pen_dirind, fwd_bwd_average, map_dpar_bwd, map_dpar_fwd) & !$omp shared(mesh, map, penalisation, values, vals_fwd, vals_bwd, gradients, & !$omp penalised_values, staggered) ! Compute map values from forward and backward plane if (staggered) then call map%upstream_stag_from_stag_fwd(values%hfwd, vals_fwd) call map%upstream_stag_from_stag_bwd(values%hbwd, vals_bwd) else call map%upstream_cano_from_cano_fwd(values%hfwd, vals_fwd) call map%upstream_cano_from_cano_bwd(values%hbwd, vals_bwd) endif !$omp do do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) if (present(gradients)) then grad = gradients(ki) else grad = 0.0_GP endif if (staggered) then map_dpar_bwd = map%dpar_stag_stag_bwd(l) map_dpar_fwd = map%dpar_stag_stag_fwd(l) else map_dpar_bwd = map%dpar_cano_cano_bwd(l) map_dpar_fwd = map%dpar_cano_cano_fwd(l) endif pen_dirind = penalisation%get_dirindfun_val(ki) fwd_bwd_average = (1.0_GP - abs(pen_dirind)) * 0.5_GP * (vals_bwd(l) + vals_fwd(l)) if (pen_dirind > 0.0_GP) then ! if magnetic field directed towards target use back map values penalised_values(ki) = + pen_dirind*(vals_bwd(l) + grad * map_dpar_bwd) & + fwd_bwd_average !Correction for the case where abs(pen_dirind) isn't 1.0 else ! if magnetic field directed away from target use forward map values penalised_values(ki) = - pen_dirind*(vals_fwd(l) - grad * map_dpar_fwd) & + fwd_bwd_average !Correction for the case where abs(pen_dirind) isn't 1.0 endif enddo !$omp end do !$omp end parallel end subroutine subroutine apply_none(mesh_any, var, penalised_values) !! Skip penalisation application type(mesh_cart_t), intent(in) :: mesh_any !! Mesh (cano or stag) type(variable_t), intent(in) :: var !! Input variable real(GP), intent(out), dimension(mesh_any%get_n_points_inner()) :: penalised_values !! Array on unaltered inner-mesh values integer :: l, ki !$omp parallel default(none) private(l, ki) & !$omp shared(mesh_any, var, penalised_values) !$omp do do ki = 1, mesh_any%get_n_points_inner() l = mesh_any%inner_indices(ki) penalised_values(ki) = var%vals(l) enddo !$omp end do !$omp end parallel end subroutine end module penalisation_neutrals_m