module penalisation_m !! Penalisation use MPI use netcdf use precision_grillix_m, only : GP, GP_NAN, GP_EPS use comm_handler_m, only : comm_handler_t use parallelisation_setup_m, only : is_rank_info_writer use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use error_handling_grillix_m, only: handle_error_netcdf, handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER, & GRILLIX_ERR_NAMELIST, GRILLIX_ERR_ALLOC, GRILLIX_ERR_SOLVER2D use descriptors_braginskii_m, only : BND_TYPE_NEUMANN use descriptors_m, only : DISTRICT_WALL, DISTRICT_OUT use boundaries_perp_m, only : set_perpbnds, extrapolate_ghost_points use mesh_index_helpers_m, only : find_neighbor_inds use parallel_map_m, only : parallel_map_t use equilibrium_storage_m, only : equilibrium_storage_t use inplane_operators_m, only : inplane_operators_t ! from PARALLAX use screen_io_m, only : get_stdout, progress_bar use constants_m, only : TWO_PI use descriptors_m, only : DISTRICT_OUT use fieldline_tracer_m, only : trace use elementary_functions_m, only : step_tanh, step_hermite, heaviside use circular_equilibrium_m, only : circular_t use root_finding_m, only : func_1d_t, find_zero use mpi_mapping_auxiliaries_m, only : getdata_fwdbwdplane use multigrid_m, only : multigrid_t use helmholtz_solver_factory_m, only : parameters_helmholtz_solver_factory, helmholtz_solver_factory use helmholtz_solver_m, only : helmholtz_solver_t implicit none integer, private, parameter :: PENFUNS_TYPE_TANH = 464 !! Descriptor for penalisation function type: tanh integer, private, parameter :: PENFUNS_TYPE_HERMITE = 465 !! Descriptor for penalisation function type: hermite polynomials integer, private, parameter :: PEN_METHOD_VIA_TRACE = 811 !! Descriptor for penalisation method: via_trace integer, private, parameter :: PEN_METHOD_VIA_IN_VESSEL = 812 !! Descriptor for penalisation method: via_in_vessel integer, private, parameter :: PEN_METHOD_NONE = 813 !! Descriptor for penalisation method: NONE integer, private, parameter :: PEN_METHOD_VIA_SUBTRACE = 814 !! Descriptor for penalisation method: via_subtrace integer, private, parameter :: PEN_METHOD_VIA_STABLE_TRACE = 815 !! Descriptor for penalisation method: via_stable_trace integer, private, parameter :: PEN_METHOD_VIA_RHO = 816 !! Descriptor for penalisation method: via_rho type, private, extends(func_1d_t) :: func_1d_charfun_t !! Test function type for setting charfun at target contains procedure :: func => func_hermite_aux end type type, public :: penalisation_t !! Datatype for penalisation integer, private :: pen_method = PEN_METHOD_VIA_TRACE !! Method for penalisation, options: !! - VIA_TRACE : Based on field line tracing with smooth transition !! - VIA_IN_VESSEL: Based on equilibrium's in_vessel function, hard transition !! - VIA_RHO: Based on the rho flux surface label integer, private :: penfuns_type !! Type of penalisation function integer, private :: hermite_order !! Order of Hermite polynomial used as penalisation function real(GP), private :: charfun_parwidth !! Parallel decay length of characteristic penalisation function real(GP), private :: charfun_radlimwidth !! Radial decay length (only for limiter geometry) real(GP), private :: dirindfun_parwidth !! Parallel decay length of direction indicator function real(GP), public, allocatable, dimension(:) :: charfun !! Characteristic function of penalisation !! Formerly known as chi real(GP), public, allocatable, dimension(:) :: dirindfun !! Function indicating direction of magnetic field (towards/away from target) !! Formerly known as xi or zeta real(GP), private :: dphi !! Toroidal grid distance real(GP), private :: dphi_max !! Maximum angle to be traced for building penalisation real(GP), private :: charfun_at_target !! Contourvalue of characteristic function at location of target plate real(GP), private :: max_step_size !! Maximum step size for the tracing real(GP), private :: rho_min !! For the build VIA_STABLE_TRACE: grid points that have smaller rho !! values than this, all penalization functions are automatically zero integer, public, allocatable, dimension(:) :: p_inds !! Indices on full mesh which lie in penalisation region integer, public :: n_p_inds !! Size of p_inds integer, public, allocatable, dimension(:) :: pb_inds !! Indices on full mesh with charfun value 0 that !! have direct neighbors in the penalisation region integer, public :: n_pb_inds !! Size of pb_inds contains procedure, public :: set_parameters => set_parameters_penalisation ! >>>>>>>>>>> ! Build routines procedure, public :: build => build_penalisation procedure, public :: build_wip3d => build_wip3d_penalisation procedure, private :: build_via_trace procedure, private :: build_via_subtrace procedure, private :: build_via_in_vessel procedure, private :: build_via_rho procedure, private :: build_pen_inds procedure, private :: build_via_stable_trace final :: destructor ! >>>>>>>>>>> ! I/O routines procedure, public :: write_netcdf => write_netcdf_penalisation procedure, public :: read_netcdf => read_netcdf_penalisation procedure, public :: display => display_penalisation ! >>>>>>>>>>> ! Will probably be deprecated procedure, public :: get_dphi procedure, public :: get_dphi_max procedure, public :: get_penfuns_type procedure, public :: get_hermite_order procedure, public :: get_charfun_parwidth procedure, public :: get_charfun_radlimwidth procedure, public :: get_dirindfun_parwidth procedure, public :: get_charfun_at_target procedure, public :: get_charfun_val procedure, public :: get_dirindfun_val end type interface real(GP) module function func_hermite_aux(this, t, iuser, ruser) !! Wrapper function for zero-finding routine class(func_1d_charfun_t), intent(in) :: this !! Instance of test function class real(GP), intent(in) :: t ! Abscissa integer, dimension(:), intent(in) :: iuser !! Integer inputs, !! (1) sets hermite_order real(GP), dimension(:), intent(in) :: ruser !! Float inputs, !! (1) sets charfun_at_target !! (2) sets charfun_parwidth end function end interface interface module subroutine set_parameters_penalisation(self, dphi, filename, & pen_method_in, & penfuns_type_in, hermite_order_in, & charfun_parwidth_in, charfun_radlimwidth_in, & dirindfun_parwidth_in, dphi_max_in, charfun_at_target_in, & max_step_size_in, rho_min_in) !! Sets parameters for penalisation, either via namelist from file, or setting explicitly implicit none class(penalisation_t), intent(inout) :: self !! Instance of the type real(GP), intent(in) :: dphi !! Toroidal grid distance character(*), intent(in), optional :: filename !! Filename where parameter are read from integer, intent(in), optional :: pen_method_in !! Method of penalisation integer, intent(in), optional :: penfuns_type_in !! Type of penalisation function integer, intent(in), optional :: hermite_order_in !! Order of Hermite polynomial used as penalisation function real(GP), intent(in), optional :: charfun_parwidth_in !! Parallel decay length of characteristic penalisation function real(GP), intent(in), optional :: charfun_radlimwidth_in !! Radial decay length (only for limiter geometry) real(GP), intent(in), optional :: dirindfun_parwidth_in !! Parallel decay length of direction indicator function real(GP), intent(in), optional :: dphi_max_in !! Parallel decay length of direction indicator function real(GP), intent(in), optional :: charfun_at_target_in !! Contourvalue of characteristic function at location of target plate real(GP), intent(in), optional :: max_step_size_in !! Maximum step size for the tracing real(GP), intent(in), optional :: rho_min_in !! For the build VIA_STABLE_TRACE: grid points that have smaller rho !! values than this, all penalization functions are automatically zero end subroutine module subroutine build_penalisation(self, comm_handler, equi, mesh, multigrid, dbgout) !! Builds penalisation functions class(penalisation_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Comunication handler class(equilibrium_t) :: equi !! Equilibrium (not changed) type(mesh_cart_t), intent(inout) :: mesh !! Mesh type(multigrid_t), intent(inout) :: multigrid !! Multigrid integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine build_via_trace(self, comm_handler, equi, mesh, dbgout) !! Builds penalisation functions based on field line tracing !! with smooth transition class(penalisation_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Comunication handler class(equilibrium_t) :: equi !! Equilibrium (not changed) type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine build_via_subtrace(self, comm_handler, equi, mesh, multigrid, dbgout) !! Builds penalisation functions based on sub-field line tracing !! subsequently from plane to plane with smooth transition class(penalisation_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Comunication handler class(equilibrium_t) :: equi !! Equilibrium (not changed) type(mesh_cart_t), intent(inout) :: mesh !! Mesh type(multigrid_t), intent(inout) :: multigrid !! Multigrid integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine build_via_in_vessel(self, comm_handler, equi, mesh, dbgout) !! Builds penalisation functions based on !! equilibrium's in_vessel function (hard transition) class(penalisation_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Comunication handler class(equilibrium_t) :: equi !! Equilibrium (not changed) type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine build_via_stable_trace(self, comm_handler, equi, mesh, dbgout) !! Builds penalisation functions based on field line tracing that can handle !! every target shape, with smooth transition class(penalisation_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Comunication handler class(equilibrium_t) :: equi !! Equilibrium (not changed) type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine build_via_rho(self, comm_handler, filepath, equi, mesh, dbgout) !! Builds penalisation functions based on !! equilibrium's rho flux surface label class(penalisation_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Comunication handler character(len=*), intent(in) :: filepath !! Path to file where to read parameters specific for rho penalisation from class(equilibrium_t) :: equi !! Equilibrium (not changed) type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine build_pen_inds(self, mesh, dbgout) !! Build storage of indices in penalisation/bordering penalisation region !! Note that index is w.r.t full mesh, but points must be from inner mesh class(penalisation_t), intent(inout) :: self !! Instance of class type(mesh_cart_t) :: mesh !! Mesh integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine display_penalisation(self) !! Displays information on penalisation class(penalisation_t), intent(inout) :: self !! Instance of the type end subroutine module subroutine write_netcdf_penalisation(self, fgid) !! Writes penalisation to netcdf file class(penalisation_t), intent(in) :: self !! Instance of the type integer, intent(in) :: fgid !! netcdf file or group id end subroutine module subroutine read_netcdf_penalisation(self, fgid) !! Reads penalisation from netcdf file class(penalisation_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: fgid !! netcdf file or group id end subroutine module subroutine destructor(self) !! Frees memory associated with penalisation type(penalisation_t), intent(inout) :: self !! Instance of the type end subroutine module subroutine build_wip3d_penalisation(pen_cano, pen_stag, comm_handler, & equi, equi_on_cano, equi_on_stag, & mesh_cano, mesh_stag, multigrid_cano, multigrid_stag, & map, dbgout) !! Builds penalisation functions for 3d equilibria based on parallel diffusion equation !! This is yet an experimental feature class(penalisation_t), intent(inout) :: pen_cano !! Penalisation for canonical mesh class(penalisation_t), intent(inout) :: pen_stag !! Penalisation for staggered mesh type(comm_handler_t), intent(in) :: comm_handler !! Comunication handler class(equilibrium_t) :: equi !! Equilibrium (not changed) class(equilibrium_storage_t), intent(inout) :: equi_on_cano !! Equilibrim quantities on canonical mesh class(equilibrium_storage_t), intent(inout) :: equi_on_stag !! Equilibrim quantities on staggered mesh type(mesh_cart_t), intent(inout) :: mesh_cano !! Mesh (canonical) type(mesh_cart_t), intent(inout) :: mesh_stag !! Mesh (staggered) type(multigrid_t), intent(inout) :: multigrid_cano !! Multigrid (canonical) type(multigrid_t), intent(inout) :: multigrid_stag !! Multigrid (staggered) type(parallel_map_t), intent(in) :: map !! Map integer, intent(in), optional :: dbgout !! Debug output level end subroutine end interface contains pure integer function get_penfuns_type(self) !! Returns type of penalisation function class(penalisation_t), intent(in) :: self !! Instance of the type get_penfuns_type = self%penfuns_type end function pure integer function get_hermite_order(self) !! Returns order of Hermite polynomial used as penalisation function class(penalisation_t), intent(in) :: self !! Instance of the type get_hermite_order = self%hermite_order end function pure real(GP) function get_charfun_parwidth(self) !! Returns parallel decay length of characteristic penalisation function class(penalisation_t), intent(in) :: self !! Instance of the type get_charfun_parwidth = self%charfun_parwidth end function pure real(GP) function get_charfun_radlimwidth(self) !! Radial radial decay length (only for limiter geometry) class(penalisation_t), intent(in) :: self !! Instance of the type get_charfun_radlimwidth = self%charfun_radlimwidth end function pure real(GP) function get_dirindfun_parwidth(self) !! Returns parallel decay length of direction indicator function class(penalisation_t), intent(in) :: self !! Instance of the type get_dirindfun_parwidth = self%dirindfun_parwidth end function pure real(GP) function get_dphi(self) !! Toroidal grid distance class(penalisation_t), intent(in) :: self !! Instance of the type get_dphi = self%dphi end function pure real(GP) function get_dphi_max(self) !! Maximum angle to be traced for building penalisation class(penalisation_t), intent(in) :: self !! Instance of the type get_dphi_max = self%dphi_max end function pure real(GP) function get_charfun_at_target(self) !! Contourvalue of characteristic function at target plate class(penalisation_t), intent(in) :: self !! Instance of the type get_charfun_at_target = self%charfun_at_target end function pure real(GP) function get_charfun_val(self, ki) !! Returns value of characteristic function class(penalisation_t), intent(in) :: self !! Instance of the type integer, intent(in) :: ki !! Running index for inner grid get_charfun_val = self%charfun(ki) end function pure real(GP) function get_dirindfun_val(self, ki) !! Returns value of direction indicator function class(penalisation_t), intent(in) :: self !! Instance of the type integer, intent(in) :: ki !! Running index for inner grid get_dirindfun_val = self%dirindfun(ki) end function end module