module equilibrium_m use precision_m, only: FP implicit none private ! Base case for equilibria ! "abstract" type since we will never make a general equilibrium ! All of the methods are "deferred", which is to say that their implementation will only be specified ! in subclasses of equilibrium type,public,abstract :: equilibrium_t !! Class-instance variables logical, public :: initialized = .false. real(kind=FP), public :: x0 !! Magnetic axis x = R/R0 (in normalised units) real(kind=FP), public :: y0 !! Magnetic axis y = Z/R0 (in normalised units) real(FP), public :: phi0 = 0.0_FP !! Magnetic axis phi real(kind=FP), public :: xmin, xmax, ymin, ymax !! Box limits real(kind=FP), public :: rhomax, rhomin !! Global limits for rho (rho = normalised psi, n.b. there may also be region-specific limits defined in equi) contains ! Non-deferred procedures are shared by all equilibrium objects procedure, public, pass(self) :: absb procedure, public, pass(self) :: bpol ! Class-instance procedures which must be declared for each sub-class of the equilibria base class ! N.b. they are "deferred" -- implementations will be given in subclasses procedure(init), public, deferred, pass(self) :: init procedure(display), public, deferred, pass(self) :: display procedure(debug), public, deferred, pass(self) :: debug procedure(is_axisymmetric), public, deferred, pass(self) :: is_axisymmetric procedure(rho), public, deferred, pass(self) :: rho procedure(bx), public, deferred, pass(self) :: bx procedure(by), public, deferred, pass(self) :: by procedure(btor), public, deferred, pass(self) :: btor procedure(jacobian), public, deferred, pass(self) :: jacobian procedure(epol), public, deferred, pass(self) :: epol procedure(erad), public, deferred, pass(self) :: erad procedure(district), public, deferred, pass(self) :: district procedure(in_vessel), public, deferred, pass(self) :: in_vessel procedure(mag_axis_loc), public, deferred, pass(self) :: mag_axis_loc end type equilibrium_t ! An "interface" for functions defines the type and intent of arguments, without specifying the routine itself ! All equilibrium objects must have the following procedures abstract interface ! N.b. due to the call to the interpolator, `self` must be inout subroutine init(self, filename) !! Initialises equilibrium, i.e.~reads any required parameters from file. import :: equilibrium_t class(equilibrium_t), intent(inout) :: self character(len = *), intent(in), optional :: filename end subroutine init subroutine display(self) !! Print to console minimal information about the equilibrium import :: equilibrium_t class(equilibrium_t), intent(in) :: self end subroutine display subroutine debug(self) !! Print to console extensive information about the equilibrium import :: equilibrium_t class(equilibrium_t), intent(in) :: self end subroutine debug logical function is_axisymmetric(self) !! Returns if equilibrium is axisymmetric or 3D (stellarator) import :: equilibrium_t class(equilibrium_t), intent(in) :: self end function real(kind=FP) function rho(self, x, y, phi) !! Flux surface label normalised such that rho = 0 at the magnetic !! axis, and rho = 1 at the separatrix import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) end function rho real(kind=FP) function bx(self, x, y, phi) !! Radial field component normalized to on axis field strength. import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) end function bx real(kind=FP) function by(self, x, y, phi) !! Vertical field component normalized to on axis field strength. import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) end function by real(kind=FP) function btor(self, x, y, phi) !! Toroidal field component normalized to on axis field strength. import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) end function btor real(kind=FP) function jacobian(self, x, y, phi) !! Jacobian of underlying coordinate system import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) end function jacobian subroutine epol(self, x, y, phi, epolx, epoly) !! Poloidal (along flux surface) unit vector. import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(kind=FP), intent(out) :: epolx, epoly end subroutine epol subroutine erad(self, x, y, phi, eradx, erady) !! Radial (across flux surface) unit vector. import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(kind=FP), intent(out) :: eradx, erady end subroutine erad integer function district(self, x, y, phi) !! Returns a descriptor which specifies whether a given point is in the closed, open or private flux region, !! and whether that point is within the computational domain import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) end function district logical function in_vessel(self, x, y, phi) !! Returns a boolean of whether a given point is within the reactor vessel (inside the bounding surface of !! the first-wall and divertor) or outside this bounding surface. !! N.b. DISTRICT_CORE is still within the first-wall/divertor import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) end function in_vessel subroutine mag_axis_loc(self, phi, axis_x, axis_y) !! Calculated the coordinates of the magnetic axis import :: FP, equilibrium_t class(equilibrium_t), intent(in) :: self !! Instance of class real(kind=FP), intent(in) :: phi !! Toroidal angle real(FP), intent(out) :: axis_x !! x-coordinate of the magnetic axis real(FP), intent(out) :: axis_y !! y-coordinate of the magnetic axis end subroutine mag_axis_loc end interface contains real(kind=FP) function absb(self, x, y, phi) !! Absolute value of magnetic field. class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) absb = sqrt(self%bx(x, y, phi)**2 + self%by(x, y, phi)**2 & + self%btor(x, y, phi)**2) end function absb real(kind=FP) function bpol(self, x, y, phi) !! Magnetic field component b poloidal normalised to absolute value of B (on axis) class(equilibrium_t), intent(in) :: self real(kind=FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) bpol = sqrt(self%bx(x, y, phi)**2 + self%by(x, y, phi)**2) end function bpol end module equilibrium_m