module inplane_operators_m !! Operators acting within poloidal plane (~perpendicular) use precision_grillix_m, only : GP, GP_NAN use screen_io_m, only : get_stdout use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use equilibrium_storage_m, only : equilibrium_storage_t use slab_equilibrium_m, only: slab_t use circular_equilibrium_m, only: circular_t use dommaschk_equilibrium_m, only: dommaschk_t implicit none type, public :: inplane_operators_t !! Datatype for inplane operators class(equilibrium_t), private, pointer :: equi => null() !! Pointer to equilibrium class(mesh_cart_t), private, pointer :: mesh => null() !! Pointer to the mesh class(equilibrium_storage_t), private, pointer :: equi_on_mesh => null() !! Pointer to the equilibrium quantities on mesh integer, private :: n_points !! Number of mesh points integer, private :: n_points_inner !! Number of inner mesh points real(GP), private :: hf_inv = GP_NAN !! Inverse of normalised grid spacing real(GP), private :: hf_inv_sqrd = GP_NAN !! Inverse of normalised grid spacing squared real(GP), private :: delta = GP_NAN !! Normalisation factor R0 / rhos for curvature operator real(GP), private :: deltabeta = GP_NAN !! Normalisation factor delta*beta for flutter operators real(GP), private, allocatable, dimension(:) :: wrk1 !! Work arrays for subroutine apply_dissipation_inplane !! Need to be defined here such that their scope is shared !! apply_dissipation_inplane can be safely and correctly called from a parallel region real(GP), private, allocatable, dimension(:) :: wrk2 !! Another work array procedure(curvature_interface), public, pointer :: curvature => null() !! Curvature operator procedure(curvature_interface), public, pointer :: curvature_sqrd => null() !! Curvature operator contains procedure, public :: init => init_inplane_operators procedure, public :: ddx procedure, public :: ddy procedure, public :: nabla_pol_sqrd procedure, public :: arakawa procedure, public :: doubled_arakawa procedure, public :: curvature_btor_based procedure, public :: curvature_ddy_based procedure, public :: curvature_sqrd_btor_based procedure, public :: curvature_sqrd_ddy_based procedure, public :: polarisation_advection procedure, public :: apply_dissipation_inplane procedure, private :: arakawa_loc procedure, private :: doubled_arakawa_loc procedure, public :: flutter_grad procedure, public :: flutter_div final :: destructor end type interface module subroutine init_inplane_operators(self, equi, mesh, equi_on_mesh, norm_scale, norm_flutter) !! Initialises inplane operators class(inplane_operators_t), intent(inout) :: self !! Instance of the type class(equilibrium_t), target, intent(inout) :: equi !! Equilibrium type(mesh_cart_t), target, intent(inout) :: mesh !! Mesh class(equilibrium_storage_t), target, intent(inout) :: equi_on_mesh !! Equilibrim quantities on mesh real(GP), intent(in) :: norm_scale ! Some normalisation scale (~rhos/R_0) real(GP), intent(in) :: norm_flutter ! Normalisation of flutter terms (~beta*delta) end subroutine pure real(GP) module function ddx(self, u, l) !! Evaluates derivative in x direction class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which x-derivative shall be computed end function pure real(GP) module function ddy(self, u, l) !! Evaluates derivative in y direction class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which y-derivative shall be computed end function real(GP) module function nabla_pol_sqrd(self, u, l, co) !! Evaluates J^(-1) d/dx(J co du/dx u) + d/dy(co du/dy) !! where J is the Jacobian class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which nabla_pol_sqrd is be computed real(GP), dimension(self%n_points), optional, intent(in) :: co !! Coefficient, if not present it is assumed 1 end function pure real(GP) module function arakawa(self, u, v, l) !! Evaluates Arakawa bracket [u,v] = du/dx*dv/dy - du/dy*dv/dx !! second order, see [A. Arakawa, J. Comput. Phys. 135, 103 (1997)] class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u real(GP), dimension(self%n_points), intent(in) :: v !! Values of variable v integer, intent(in) :: l !! Index, for which Arakawa Bracket is be computed end function pure real(GP) module function doubled_arakawa(self, u, v, l) !! Evaluates doubled Arakawa bracket [u,[u,v]], second order class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u real(GP), dimension(self%n_points), intent(in) :: v !! Values of variable v integer, intent(in) :: l !! Index, for which double Arakawa bracket is be computed end function real(GP) module function curvature_interface(self, u, l) !! Interface for (squared) curvature operator class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which curvature is computed end function real(GP) module function curvature_btor_based(self, u, l) !! Curvature operator !! Computation based on toroidal magnetic field class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which curvature is computed end function pure real(GP) module function curvature_ddy_based(self, u, l) !! Curvature operator !! Simple computation via vertical derivative -2 d/dy !! Useful for slab/circular geometry class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which curvature is computed end function real(GP) module function curvature_sqrd_btor_based(self, u, l) !! Evaluates squared curvature operator !! Computation based on toroidal magnetic field class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which curvature is computed end function pure real(GP) module function curvature_sqrd_ddy_based(self, u, l) !! Evaluates squared curvature operator !! Simple computation via vertical derivative 4 d^2/dy^2 !! second order, central finite difference class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which curvature is computed end function real(GP) module function polarisation_advection(self, dia_fac, co, adv_pot, potv, tiv, nev, l) !! Evaluates non-linear advection div(co [adv_pot, v]_{x,y}) !! where v = grad(pot) + diafac*grad(pi) / ne class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), intent(in) :: dia_fac !! Factor in front of diamagnetic term real(GP), dimension(self%n_points), intent(in) :: co !! Values of coefficient real(GP), dimension(self%n_points), intent(in) :: adv_pot !! Values of the first input in bracket [adv_pot, v] real(GP), dimension(self%n_points), intent(in) :: potv !! Values of electrostatic potential real(GP), dimension(self%n_points), intent(in) :: tiv !! Values of ion temperature real(GP), dimension(self%n_points), intent(in) :: nev !! Values of density integer, intent(in) :: l !! Index, for which term is computed end function module subroutine apply_dissipation_inplane(self, u, ndiss, du, co) !! Applies dissipation of order ndiss: !! du = div (co grad (nabla_pol_sqrd^{ndiss-1} u) ) class(inplane_operators_t), intent(inout) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: ndiss !! Order of dissipation operator real(GP), dimension(self%n_points_inner), intent(out) :: du !! Change added with dissipation real(GP), dimension(self%n_points), optional, intent(in) :: co !! Coefficient applied in last Laplacian operation end subroutine pure real(GP) module function arakawa_loc(self, zeta, psi) !! Evaluates arakawa bracket on local 3x3 stencil without its prefactor 1/(12 hf^2) class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(-1:1,-1:1), intent(in) :: zeta !! First entry in arakawa bracket given as local 3x3 stencil real(GP), dimension(-1:1,-1:1), intent(in) :: psi !! First entry in arakawa bracket given as local 3x3 stencil end function pure real(GP) module function doubled_arakawa_loc(self, u_loc, v_loc, jac_loc) !! Evaluates arakawa bracket on local 5x5 stencil !! If jac_loc not present: [u,[u,v]] !! If jac_loc present: 1/J*[J*u, 1/J*[J*u,v]] class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(-2:2,-2:2), intent(in) :: u_loc !! Values of variable u given as local 5x5 stencil real(GP), dimension(-2:2,-2:2), intent(in) :: v_loc !! Values of variable v given as local 5x5 stencil real(GP), dimension(-2:2,-2:2), intent(in), optional :: jac_loc !! Jacobian given as local 5x5 stencil end function pure real(GP) module function flutter_grad(self, apar, u, l) !! Evaluates flutter gradient of f !! using Arakawa bracket co*[apar,u]/btor class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: apar !! Values of variable apar real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which flutter_grad is be computed end function pure real(GP) module function flutter_div(self, apar, u, l) !! Evaluates flutter divergence of f !! using Arakawa bracket co*[apar,u]/btor class(inplane_operators_t), intent(in) :: self !! Instance of the type real(GP), dimension(self%n_points), intent(in) :: apar !! Values of variable apar real(GP), dimension(self%n_points), intent(in) :: u !! Values of variable u integer, intent(in) :: l !! Index, for which flutter_div is be computed end function module subroutine destructor(self) !! Frees memory associated with karniadakis type(inplane_operators_t), intent(inout) :: self !! Instance of the type end subroutine end interface end module