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       !< boundary conditions
    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, allocatable    :: masses(:)        !< 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, allocatable         :: n_ghost(:)   !< 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, private     :: dim           !< spatial dimension
    integer              :: center        !< index to the central point of the stencil
    integer              :: size          !< size (number of points)
    integer              :: npoly         !< this is modified in derivatives_oct_m::derivatives_build()
    integer, allocatable :: points(:, :)  !< list of points in spatial index notation (1:dim, 1:size)

    ! 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