Derivatives
Octopus use finite differences to evaluate derivatives, such as the Laplacian (kinetic energy). The derivative at a point of the mesh is a weighted sum over neighboring points. For instance, the general form for the Laplacian is: $$ \nabla^2f(n_xh,,n_yh) = \sum_i^n\sum_j^n\frac{c_{ij}}h,f(n_xh + ih,,n_yh+jh) $$ The coefficients $(c_{ij})$ depend on the mesh and number of points used and define the stencil.
The derivatives object contains operators for the gradient and the Laplacian, as well as information about the order and the stencil type, and information about boundaries:
type derivatives_t
! Components are public by default
type(boundaries_t) :: boundaries
type(mesh_t), pointer :: mesh => NULL() !< pointer to the underlying mesh
integer :: dim = 0 !< dimensionality of the space (space%dim)
integer :: order = 0 !< order of the discretization (value depends on stencil)
integer :: stencil_type = 0 !< type of discretization
FLOAT :: masses(MAX_DIM) = M_ZERO !< we can have different weights (masses) per space direction
!> If the so-called variational discretization is used, this controls a
!! possible filter on the Laplacian.
FLOAT, private :: lapl_cutoff = M_ZERO
type(nl_operator_t), allocatable, private :: op(:) !< op(1:conf%dim) => gradient
!! op(conf%dim+1) => Laplacian
type(nl_operator_t), pointer :: lapl => NULL() !< these are just shortcuts for op
type(nl_operator_t), pointer :: grad(:) => NULL()
integer :: n_ghost(MAX_DIM) = 0 !< ghost points to add in each dimension
#if defined(HAVE_MPI)
integer, private :: comm_method = 0
#endif
type(derivatives_t), pointer :: finer => NULL()
type(derivatives_t), pointer :: coarser => NULL()
type(transfer_table_t), pointer :: to_finer => NULL()
type(transfer_table_t), pointer :: to_coarser => NULL()
end type derivatives_t
The derivative operators themselves are represented as non-local operators:
type nl_operator_t
private
type(stencil_t), public :: stencil
type(mesh_t), pointer :: mesh => NULL() !< pointer to the underlying mesh
integer, allocatable :: nn(:) !< the size of the stencil at each point (for curvilinear coordinates)
integer, public :: np = 0 !< number of points in mesh
!> When running in parallel mode, the next three arrays are unique on each node.
integer, allocatable, public :: index(:,:) !< index of the points. Unique on each parallel process.
FLOAT, allocatable, public :: w(:,:) !< weights. Unique on each parallel process.
logical, public :: const_w = .true. !< are the weights independent of index i
type(accel_mem_t), public :: buff_weights !< buffer with constant weights
type(accel_mem_t), public :: buff_half_weights !< buffer with weights multiplied by -1/2
character(len=40) :: label
!> the compressed index of grid points
integer, public :: nri = 0
integer, allocatable, public :: ri(:,:)
integer, allocatable, public :: rimap(:)
integer, allocatable, public :: rimap_inv(:)
integer :: ninner = 0
integer :: nouter = 0
type(nl_operator_index_t) :: inner
type(nl_operator_index_t) :: outer
type(accel_kernel_t) :: kernel
type(accel_mem_t) :: buff_imin
type(accel_mem_t) :: buff_imax
type(accel_mem_t) :: buff_ri
type(accel_mem_t) :: buff_map
type(accel_mem_t) :: buff_all
type(accel_mem_t) :: buff_inner
type(accel_mem_t) :: buff_outer
type(accel_mem_t) :: buff_stencil
type(accel_mem_t) :: buff_ip_to_xyz
type(accel_mem_t) :: buff_xyz_to_ip
end type nl_operator_t
type stencil_t
! Components are public by default
integer :: center
integer :: size
integer :: npoly
integer, allocatable :: points(:, :)
! The stargeneral arms
type(stargeneral_arms_t) :: stargeneral
end type stencil_t
type stargeneral_arms_t
! Components are public by default
integer :: arms(1:3,1:3) = reshape([0,0,0,0,0,0,0,0,0], (/3, 3/))
integer :: narms = 0
end type stargeneral_arms_t