equilibrium_storage_m.f90 Source File

module equilibrium_storage_m
    !! Storage of equilibrium data on mesh enabling lookup instead of computation 
    !! Can be used for faster performance at certain locations 
    use netcdf
    use error_handling_grillix_m, only: handle_error_netcdf, handle_error
    use precision_grillix_m, only : GP
    use screen_io_m, only :  get_stdout
    use equilibrium_m, only : equilibrium_t
    use mesh_cart_m, only : mesh_cart_t
    implicit none
    
        
    type, public ::  equilibrium_storage_t
        !! Datatype for storage of equilibrium on mesh
        real(GP), allocatable, dimension(:), public :: rho
        !! rho on mesh
        real(GP), allocatable, dimension(:), public :: bx
        !! bx on mesh
        real(GP), allocatable, dimension(:), public :: by
        !! by on mesh
        real(GP), allocatable, dimension(:), public :: btor
        !! btor on mesh
        real(GP), allocatable, dimension(:,:), public :: epol
        !! epol on mesh
        real(GP), allocatable, dimension(:,:), public :: erad
        !! erad on mesh
        integer, allocatable, dimension(:), public :: district
        !! District on mesh
    contains
        procedure, public :: absb
        procedure, public :: fill_storage
        procedure, public :: write_netcdf => write_netcdf_equilibrium_storage
        procedure, public :: read_netcdf => read_netcdf_equilibrium_storage
        final :: destructor    
    end type

    interface
        module subroutine write_netcdf_equilibrium_storage(self, fgid)
            !! Writes stored quantities of grid to netcdf file
            class(equilibrium_storage_t), intent(in) :: self
            !! Instance of class
            integer, intent(in) :: fgid
            !! File or group id
        end subroutine

        module subroutine read_netcdf_equilibrium_storage(self, fgid)
            !! Writes stored quantities of grid to netcdf file
            class(equilibrium_storage_t), intent(inout) :: self
            !! Instance of class
            integer, intent(in) :: fgid
            !! File or group id
        end subroutine
    end interface

contains

    pure real(GP) function absb(self, l)
        !! Computes absolute value of magnetic field from stured quantities
        class(equilibrium_storage_t), intent(in) :: self
        !! Instance of class
        integer, intent(in) :: l
        !! Index
        absb = sqrt(self%bx(l)**2 + self%by(l)**2 + self%btor(l)**2)
    end function
    

    subroutine fill_storage(self, equi, mesh)
        !! Fills storage with values on mesh points
        class(equilibrium_storage_t), intent(inout) :: self
        !! Instance of class
        class(equilibrium_t), intent(inout), target :: equi
        !! Equilibrium
        type(mesh_cart_t), intent(in) :: mesh
        !! Mesh
        
        integer :: np, l
        real(GP) :: x, y, phi
                
        np = mesh%get_n_points()
        
        allocate(self%rho(np))
        allocate(self%bx(np))
        allocate(self%by(np))
        allocate(self%btor(np))
        allocate(self%epol(np,2))
        allocate(self%erad(np,2))
        allocate(self%district(np))
        
        phi = mesh%get_phi()
        !$OMP PARALLEL PRIVATE(l, x, y)
        !$OMP DO
        do l = 1, mesh%get_n_points()
            x = mesh%get_x(l) 
            y = mesh%get_y(l)
            self%rho(l) = equi%rho(x, y, phi)
            if (equi%is_axisymmetric()) then            
                call equi%epol(x, y, phi, self%epol(l,1), self%epol(l,2))
                call equi%erad(x, y, phi, self%erad(l,1), self%erad(l,2))
            else
                self%epol(l,:) = 0.0_GP ! TODO: DOMMASCHK
                self%erad(l,:) = 0.0_GP ! TODO: DOMMASCHK
            endif
            self%bx(l)    = equi%bx(x, y, phi)
            self%by(l)    = equi%by(x, y, phi)
            self%btor(l)  = equi%btor(x, y, phi)

            self%district(l) = equi%district(x, y, phi)
        enddo
        !$OMP END DO
        !$OMP END PARALLEL
        
    end subroutine
    
    subroutine destructor(self)
        !! Destructor
        type(equilibrium_storage_t), intent(inout) :: self
    
        if (allocated(self%rho))  deallocate(self%rho)
        if (allocated(self%bx))   deallocate(self%bx)
        if (allocated(self%by))   deallocate(self%by)
        if (allocated(self%btor)) deallocate(self%btor)
        if (allocated(self%epol)) deallocate(self%epol)
        if (allocated(self%erad)) deallocate(self%erad)
        if (allocated(self%district)) deallocate(self%district)

    end subroutine
    
end module