derivatives.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
 
module derivatives_oct_m
  use accel_oct_m
  use affine_coordinates_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use boundaries_oct_m
  use debug_oct_m
  use cartesian_oct_m
  use coordinate_system_oct_m
  use curv_briggs_oct_m
  use curv_gygi_oct_m
  use curv_modine_oct_m
  use global_oct_m
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use lalg_adv_oct_m
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use nl_operator_oct_m
  use par_vec_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use stencil_oct_m
  use stencil_cube_oct_m
  use stencil_star_oct_m
  use stencil_starplus_oct_m
  use stencil_stargeneral_oct_m
  use stencil_variational_oct_m
  use transfer_table_oct_m
  use types_oct_m
  use utils_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                             &
    derivatives_t,                      &
    derivatives_init,                   &
    derivatives_end,                    &
    derivatives_build,                  &
    derivatives_handle_batch_t,         &
    dderivatives_batch_start,           &
    zderivatives_batch_start,           &
    dderivatives_batch_finish,          &
    zderivatives_batch_finish,          &
    dderivatives_batch_perform,         &
    zderivatives_batch_perform,         &
    dderivatives_perform,               &
    zderivatives_perform,               &
    dderivatives_lapl,                  &
    zderivatives_lapl,                  &
    derivatives_lapl_diag,              &
    dderivatives_grad,                  &
    zderivatives_grad,                  &
    dderivatives_batch_grad,            &
    zderivatives_batch_grad,            &
    dderivatives_batch_div,             &
    zderivatives_batch_div,             &
    dderivatives_div,                   &
    zderivatives_div,                   &
    dderivatives_curl,                  &
    zderivatives_curl,                  &
    dderivatives_batch_curl,            &
    zderivatives_batch_curl,            &
    dderivatives_batch_curl_from_gradient, &
    zderivatives_batch_curl_from_gradient, &
    dderivatives_partial,               &
    zderivatives_partial,               &
    derivatives_get_lapl,               &
    derivatives_get_inner_boundary_mask, &
    derivatives_lapl_get_max_eigenvalue, &
    derivates_set_coordinates_system
 
 
  integer, parameter ::     &
    DER_BC_ZERO_F    = 0,   &  !< function is zero at the boundaries
    der_bc_zero_df   = 1,   &  
    der_bc_period    = 2       
 
  integer, parameter, public ::     &
    DER_STAR         = 1,   &
    der_variational  = 2,   &
    der_cube         = 3,   &
    der_starplus     = 4,   &
    der_stargeneral  = 5
 
  integer, parameter ::     &
    BLOCKING = 1,           &
    non_blocking = 2
 
  type derivatives_t
    ! Components are public by default
    type(boundaries_t)    :: boundaries
    type(mesh_t), pointer :: mesh => null()   
    integer               :: dim = 0
    integer               :: periodic_dim = 0
    integer               :: order = 0
    integer               :: stencil_type = 0
 
    
    logical, private      :: stencil_primitive_coordinates
 
    class(coordinate_system_t), private, pointer    :: coord_system
    logical               :: remove_zero_weight_points
 
    real(real64), allocatable    :: masses(:)
 
    real(real64), private :: lapl_cutoff = m_zero
 
    type(nl_operator_t), allocatable, private :: op(:)
    type(nl_operator_t), pointer :: lapl => null()      
    type(nl_operator_t), pointer :: grad(:) => null()
 
    integer, allocatable         :: n_ghost(:)
!#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
 
  type derivatives_handle_batch_t
    private
!#ifdef HAVE_MPI
    type(par_vec_handle_batch_t) :: pv_h
!#endif
    type(derivatives_t), pointer :: der
    type(nl_operator_t), pointer :: op
    type(batch_t),       pointer :: ff
    type(batch_t),       pointer :: opff
    logical                      :: ghost_update
    logical                      :: factor_present
    real(real64)                 :: factor
  end type derivatives_handle_batch_t
 
  type(accel_kernel_t) :: kernel_uvw_xyz, kernel_dcurl, kernel_zcurl
 
contains
 
  ! ---------------------------------------------------------
  subroutine derivatives_init(der, namespace, space, coord_system, order)
    type(derivatives_t),        target, intent(inout) :: der
    type(namespace_t),                  intent(in)    :: namespace
    class(space_t),                     intent(in)    :: space
    class(coordinate_system_t), target, intent(in)    :: coord_system
    integer, optional,                  intent(in)    :: order
 
    integer :: idir
    integer :: default_stencil
    logical :: stencil_primitive_coordinates
    character(len=128) :: kernel_name
 
    push_sub(derivatives_init)
 
    ! Non-orthogonal curvilinear coordinates are currently not implemented
    assert(.not. coord_system%local_basis .or. coord_system%orthogonal)
 
    ! copy this value to my structure
    der%dim = space%dim
    der%periodic_dim = space%periodic_dim
    der%coord_system => coord_system
 
    !%Variable DerivativesStencil
    !%Type integer
    !%Default stencil_star
    !%Section Mesh::Derivatives
    !%Description
    !% Decides what kind of stencil is used, <i>i.e.</i> which points, around
    !% each point in the mesh, are the neighboring points used in the
    !% expression of the differential operator.
    !%
    !% If curvilinear coordinates are to be used, then only the <tt>stencil_starplus</tt>
    !% or the <tt>stencil_cube</tt> may be used. We only recommend the <tt>stencil_starplus</tt>,
    !% since the cube typically needs far too much memory.
    !%Option stencil_star 1
    !% A star around each point (<i>i.e.</i>, only points on the axis).
    !%Option stencil_variational 2
    !% Same as the star, but with coefficients built in a different way.
    !%Option stencil_cube 3
    !% A cube of points around each point.
    !%Option stencil_starplus 4
    !% The star, plus a number of off-axis points.
    !%Option stencil_stargeneral 5
    !% The general star. Default for non-orthogonal grids.
    !%End
    default_stencil = der_star
    if (coord_system%local_basis) default_stencil = der_starplus
    if (.not. coord_system%orthogonal) default_stencil = der_stargeneral
 
    call parse_variable(namespace, 'DerivativesStencil', default_stencil, der%stencil_type)
 
    if (.not. varinfo_valid_option('DerivativesStencil', der%stencil_type)) then
      call messages_input_error(namespace, 'DerivativesStencil')
    end if
    call messages_print_var_option("DerivativesStencil", der%stencil_type, namespace=namespace)
 
    if (coord_system%local_basis .and.  der%stencil_type < der_cube) call messages_input_error(namespace, 'DerivativesStencil')
    ! The star stencils have no off-axis points, so they cannot represent the cross derivatives
    ! of the Laplacian in non-orthogonal coordinates.
    if (.not. coord_system%orthogonal .and. &
      (der%stencil_type == der_star .or. der%stencil_type == der_variational)) then
      call messages_input_error(namespace, 'DerivativesStencil', &
        'stencil_star and stencil_variational cannot be used for non-orthogonal cells. Use stencil_stargeneral.')
    end if
    if (der%stencil_type == der_variational) then
      !%Variable DerivativesLaplacianFilter
      !%Type float
      !%Default 1.0
      !%Section Mesh::Derivatives
      !%Description
      !% Undocumented
      !%End
      call parse_variable(namespace, 'DerivativesLaplacianFilter', m_one, der%lapl_cutoff)
    end if
 
    !%Variable DerivativesOrder
    !%Type integer
    !%Default 4
    !%Section Mesh::Derivatives
    !%Description
    !% This variable gives the discretization order for the approximation of
    !% the differential operators. This means, basically, that
    !% <tt>DerivativesOrder</tt> points are used in each positive/negative
    !% spatial direction, <i>e.g.</i> <tt>DerivativesOrder = 1</tt> would give
    !% the well-known three-point formula in 1D.
    !% The number of points actually used for the Laplacian
    !% depends on the stencil used. Let <math>O</math> = <tt>DerivativesOrder</tt>, and <math>d</math> = <tt>Dimensions</tt>.
    !% <ul>
    !% <li> <tt>stencil_star</tt>: <math>2 O d + 1</math>
    !% <li> <tt>stencil_cube</tt>: <math>(2 O + 1)^d</math>
    !% <li> <tt>stencil_starplus</tt>: <math>2 O d + 1 + n</math> with <i>n</i> being 8
    !% in 2D and 24 in 3D.
    !% </ul>
    !%End
    call parse_variable(namespace, 'DerivativesOrder', 4, der%order)
    ! overwrite order if given as argument
    if (present(order)) then
      der%order = order
    end if
 
    !%Variable DerivativesRemoveZeroWeightPoints
    !%Type logical
    !%Default yes
    !%Section Mesh::Derivatives
    !%Description
    !% By default, Octopus removes the zero-weight points in the stencils.
    !% This debug variable gives the possibility to deactive this feature.
    !%End
    call parse_variable(namespace, 'DerivativesRemoveZeroWeightPoints', .true., der%remove_zero_weight_points)
 
 
#ifdef HAVE_MPI
    !%Variable ParallelizationOfDerivatives
    !%Type integer
    !%Default non_blocking
    !%Section Execution::Parallelization
    !%Description
    !% This option selects how the communication of mesh boundaries is performed.
    !%Option blocking 1
    !% Blocking communication.
    !%Option non_blocking 2
    !% Communication is based on non-blocking point-to-point communication.
    !%End
 
    call parse_variable(namespace, 'ParallelizationOfDerivatives', non_blocking, der%comm_method)
 
    if (.not. varinfo_valid_option('ParallelizationOfDerivatives', der%comm_method)) then
      call messages_input_error(namespace, 'ParallelizationOfDerivatives')
    end if
 
    call messages_obsolete_variable(namespace, 'OverlapDerivatives', 'ParallelizationOfDerivatives')
#endif
 
    !%Variable StencilPrimitiveCoordinates
    !%Type logical
    !%Section Mesh::Derivatives
    !%Description
    !% This variable controls the method for generating the stencil weights for the Laplacian.
    !% For the gradient, the weights are always computed for primitive coordinates.
    !% If set to yes, primitive coordinates are used for the polynomials and the right-hand side
    !% of the Laplacian is computed using the metric tensor and the trace of the Hessian.
    !% If set to no, Cartesian coordinates are used for the polynomials and the right-hand side
    !% of the Laplacian reduces to the Cartesian case.
    !% For some non-orthogonal grids, using primitive coordinates is necessary because the polynomials
    !% become linearly dependent, thus the corresponding matrix cannot be inverted. For some curvilinear
    !% coordinate systems, using Cartesian coordinates is more accurate.
    !%
    !% By default, use primitve coordinates except for curvilinear meshes.
    !%End
    stencil_primitive_coordinates = .not. coord_system%local_basis
    call parse_variable(namespace, 'StencilPrimitiveCoordinates', stencil_primitive_coordinates, der%stencil_primitive_coordinates)
 
    ! if needed, der%masses should be initialized in modelmb_particles_init
    safe_allocate(der%masses(1:space%dim))
    der%masses = m_one
 
    ! construct lapl and grad structures
    safe_allocate(der%op(1:der%dim + 1))
    der%grad => der%op
    der%lapl => der%op(der%dim + 1)
 
    call derivatives_get_stencil_lapl(der, der%lapl, space, coord_system)
    call derivatives_get_stencil_grad(der)
 
    ! find out how many ghost points we need in each dimension
    safe_allocate(der%n_ghost(1:der%dim))
    der%n_ghost(:) = 0
    do idir = 1, der%dim
      der%n_ghost(idir) = maxval(abs(der%lapl%stencil%points(idir, :)))
    end do
 
    nullify(der%coarser)
    nullify(der%finer)
    nullify(der%to_coarser)
    nullify(der%to_finer)
 
    if (accel_is_enabled()) then
      ! Check if we can build the uvw_to_xyz kernel
      select type (coord_system)
      type is (cartesian_t)
        ! In this case one does not need to call the kernel, so all is fine
      class default
        if (der%dim > 3) then
          message(1) = "Calculation of derivatives on the GPU with dimension > 3 are only implemented for Cartesian coordinates."
          call messages_fatal(1, namespace=namespace)
        else
          ! Each dimension adds a (source, ld, dest, ld) argument group
          write(kernel_name, '(a,i1,a)') 'uvw_to_xyz<double, Dim::D', der%dim, &
            repeat(', double*, int, double*, int', der%dim)//'>'
          call accel_kernel_start_call(kernel_uvw_xyz, 'uvw_to_xyz.cu', trim(kernel_name))
        end if
      end select
      call accel_kernel_build(kernel_dcurl, 'curl.cu', 'curl<double>')
      call accel_kernel_build(kernel_zcurl, 'curl.cu', 'curl<complex<double>>')
    end if
 
    pop_sub(derivatives_init)
  end subroutine derivatives_init
 
 
  ! ---------------------------------------------------------
  subroutine derivatives_end(der)
    type(derivatives_t), intent(inout) :: der
 
    integer :: idim
 
    push_sub(derivatives_end)
 
    assert(allocated(der%op))
 
    do idim = 1, der%dim+1
      call nl_operator_end(der%op(idim))
    end do
 
    safe_deallocate_a(der%masses)
 
    safe_deallocate_a(der%n_ghost)
 
    safe_deallocate_a(der%op)
    nullify(der%lapl, der%grad)
 
    nullify(der%coarser)
    nullify(der%finer)
    nullify(der%to_coarser)
    nullify(der%to_finer)
 
    call boundaries_end(der%boundaries)
 
    pop_sub(derivatives_end)
  end subroutine derivatives_end
 
 
  ! ---------------------------------------------------------
  subroutine derivatives_get_stencil_lapl(der, lapl, space, coord_system, name, order)
    type(derivatives_t),         intent(in)    :: der
    type(nl_operator_t),         intent(inout) :: lapl
    class(space_t),              intent(in)    :: space
    class(coordinate_system_t),  intent(in)    :: coord_system
    character(len=80), optional, intent(in)    :: name
    integer,           optional, intent(in)    :: order
 
    character(len=80) :: name_
    integer :: order_
 
    push_sub(derivatives_get_stencil_lapl)
 
    name_ = optional_default(name, "Laplacian")
    order_ = optional_default(order, der%order)
 
    ! initialize nl operator
    call nl_operator_init(lapl, name_, symm=op_symmetric)
 
    ! create stencil
    select case (der%stencil_type)
    case (der_star, der_variational)
      call stencil_star_get_lapl(lapl%stencil, der%dim, order_)
    case (der_cube)
      call stencil_cube_get_lapl(lapl%stencil, der%dim, order_)
    case (der_starplus)
      call stencil_starplus_get_lapl(lapl%stencil, der%dim, order_)
    case (der_stargeneral)
      call stencil_stargeneral_get_arms(lapl%stencil, space%dim, coord_system)
      call stencil_stargeneral_get_lapl(lapl%stencil, der%dim, order_)
    end select
 
    pop_sub(derivatives_get_stencil_lapl)
  end subroutine derivatives_get_stencil_lapl
 
 
  ! ---------------------------------------------------------
  subroutine derivatives_lapl_diag(der, lapl)
    type(derivatives_t), intent(in) :: der
    real(real64), contiguous, intent(out) :: lapl(:)
 
    push_sub(derivatives_lapl_diag)
 
    assert(ubound(lapl, dim=1) >= der%mesh%np)
 
    ! the Laplacian is a real operator
    call dnl_operator_operate_diag(der%lapl, lapl)
 
    pop_sub(derivatives_lapl_diag)
 
  end subroutine derivatives_lapl_diag
 
 
  ! ---------------------------------------------------------
  subroutine derivatives_get_stencil_grad(der)
    type(derivatives_t), intent(inout) :: der
 
 
    integer  :: ii
    character :: dir_label
 
    push_sub(derivatives_get_stencil_grad)
 
    assert(associated(der%grad))
 
    ! initialize nl operator
    do ii = 1, der%dim
      dir_label = ' '
      if (ii < 5) dir_label = index2axis(ii)
 
      call nl_operator_init(der%grad(ii), "Gradient "//dir_label, symm=op_antisymmetric)
 
      ! create stencil
      select case (der%stencil_type)
      case (der_star, der_variational)
        call stencil_star_get_grad(der%grad(ii)%stencil, der%dim, ii, der%order)
      case (der_cube)
        call stencil_cube_get_grad(der%grad(ii)%stencil, der%dim, der%order)
      case (der_starplus)
        call stencil_starplus_get_grad(der%grad(ii)%stencil, der%dim, ii, der%order)
      case (der_stargeneral)
        ! use the simple star stencil
        call stencil_star_get_grad(der%grad(ii)%stencil, der%dim, ii, der%order)
      end select
    end do
 
    pop_sub(derivatives_get_stencil_grad)
 
  end subroutine derivatives_get_stencil_grad
 
  ! ---------------------------------------------------------
  subroutine derivatives_build(der, namespace, space, mesh, qvector, regenerate, verbose)
    type(derivatives_t), intent(inout) :: der
    type(namespace_t), intent(in) :: namespace
    class(space_t), intent(in) :: space
    class(mesh_t), target, intent(in) :: mesh
    real(real64), contiguous, optional, intent(in) :: qvector(:)
    logical, optional, intent(in) :: regenerate
    logical, optional, intent(in) :: verbose
 
    integer :: i
    logical :: const_w_
    integer :: np_zero_bc
 
    push_sub(derivatives_build)
 
    if (.not. optional_default(regenerate, .false.)) then
      call boundaries_init(der%boundaries, namespace, space, mesh, qvector)
    end if
 
    assert(allocated(der%op))
    assert(der%stencil_type >= der_star .and. der%stencil_type <= der_stargeneral)
    assert(.not.(der%stencil_type == der_variational .and. mesh%use_curvilinear))
 
    der%mesh => mesh    ! make a pointer to the underlying mesh
 
    const_w_ = .true.
 
    ! need non-constant weights for curvilinear and scattering meshes
    if (mesh%use_curvilinear) const_w_ = .false.
 
    np_zero_bc = 0
 
    if (optional_default(regenerate, .false.)) then
      do i = 1, der%dim+1
        call nl_operator_end(der%op(i))
        safe_deallocate_a(der%op(i)%w)
      end do
      call derivatives_get_stencil_lapl(der, der%lapl, space, der%coord_system)
      call derivatives_get_stencil_grad(der)
    end if
 
    ! build operators
    do i = 1, der%dim+1
      call nl_operator_build(space, mesh, der%op(i), der%mesh%np, const_w = const_w_, &
        regenerate=regenerate)
      np_zero_bc = max(np_zero_bc, nl_operator_np_zero_bc(der%op(i)))
    end do
 
    assert(np_zero_bc > mesh%np .and. np_zero_bc <= mesh%np_part)
 
    select case (der%stencil_type)
 
    case (der_star, der_starplus, der_stargeneral) ! Laplacian and gradient have different stencils
      do i = 1, der%dim + 1
 
        call derivatives_make_discretization(namespace, der, 1, i, der%op(i:i), space, verbose=verbose)
      end do
 
    case (der_cube)
      ! gradients have the same stencils, so use one call to derivatives_make_discretization
      ! to solve the linear equation once for several right-hand sides
      call derivatives_make_discretization(namespace, der, der%dim, der%dim, der%op(:), space, verbose=verbose)
      ! the polynomials are evaluated differently for the Laplacian
      call derivatives_make_discretization(namespace, der, 1, der%dim+1, der%op(der%dim+1:der%dim+1), space, &
        verbose=verbose)
 
    case (der_variational)
      ! the gradients use the star stencil
      do i = 1, der%dim
        call derivatives_make_discretization(namespace, der, 1, i, der%op(i:i), space, verbose=verbose)
      end do
      ! we have the explicit coefficients for the Laplacian
      call stencil_variational_coeff_lapl(der%dim, der%order, mesh%spacing, der%lapl, alpha = der%lapl_cutoff)
      ! finalize the Laplacian as done in derivatives_make_discretization for the other stencils
      if (der%lapl%const_w) call nl_operator_build_symmetric_weights(der%lapl)
      call nl_operator_allocate_gpu_buffers(der%lapl)
      call nl_operator_update_gpu_buffers(der%lapl)
    end select
 
    pop_sub(derivatives_build)
 
  end subroutine derivatives_build
 
  ! ---------------------------------------------------------
  subroutine stencil_pol_grad(stencil_type, dim, direction, order, polynomials)
    integer, intent(in) :: stencil_type
    integer, intent(in) :: dim
    integer, intent(in) :: direction
    integer, intent(in) :: order
    integer, contiguous, intent(inout) :: polynomials(:, :)
 
    select case (stencil_type)
    case (der_star, der_stargeneral, der_variational)
      call stencil_star_polynomials_grad(direction, order, polynomials)
    case (der_starplus)
      call stencil_starplus_pol_grad(dim, direction, order, polynomials)
    case (der_cube)
      ! same stencil for gradient and Laplacian in this case
      call stencil_cube_polynomials_lapl(dim, order, polynomials)
    end select
  end subroutine stencil_pol_grad
 
  ! ---------------------------------------------------------
  subroutine stencil_pol_lapl(stencil_type, stencil, dim, order, polynomials)
    integer, intent(in) :: stencil_type
    type(stencil_t), intent(in) :: stencil
    integer, intent(in) :: dim
    integer, intent(in) :: order
    integer, contiguous, intent(inout) :: polynomials(:, :)
 
    select case (stencil_type)
    case (der_star, der_variational)
      ! the variational stencil uses the points of the star; a generic Laplacian (e.g. with a
      ! different order) is built with the star finite-difference weights
      call stencil_star_polynomials_lapl(dim, order, polynomials)
    case (der_starplus)
      call stencil_starplus_pol_lapl(dim, order, polynomials)
    case (der_stargeneral)
      call stencil_stargeneral_pol_lapl(stencil, dim, order, polynomials)
    case (der_cube)
      call stencil_cube_polynomials_lapl(dim, order, polynomials)
    end select
  end subroutine stencil_pol_lapl
 
  ! ---------------------------------------------------------
  subroutine get_rhs_lapl(coord_system, polynomials, chi, rhs, stencil_primitive_coordinates)
    class(coordinate_system_t), intent(in) :: coord_system
    integer, contiguous, intent(in) :: polynomials(:,:)
    real(real64), contiguous, intent(in) :: chi(:)
    real(real64), contiguous, intent(out) :: rhs(:)
    logical, intent(in) :: stencil_primitive_coordinates
 
    integer :: i, j, k, dim
    real(real64) :: metric_inverse(1:size(polynomials, dim=1), 1:size(polynomials, dim=1))
    real(real64) :: hessian_trace(1:size(polynomials, dim=1))
    integer :: powers(0:2)
 
    push_sub(get_rhs_lapl)
 
    dim = size(polynomials, dim=1)
 
    if (stencil_primitive_coordinates) then
      ! inverse metric and trace of Hessian in primitive coordinates
      metric_inverse = coord_system%metric_inverse(chi)
      hessian_trace = coord_system%trace_hessian(chi)
    else
      ! inverse metric and trace of Hessian in Cartesian coordinates
      metric_inverse = m_zero
      do i = 1, dim
        metric_inverse(i, i) = m_one
      end do
      hessian_trace = m_zero
    end if
 
    ! find right-hand side for operator
    rhs(:) = m_zero
    do j = 1, size(polynomials, dim=2)
      ! count the powers of the polynomials
      powers = 0
      do i = 1, dim
        if (polynomials(i, j) <= 2) then
          powers(polynomials(i, j)) = powers(polynomials(i, j)) + 1
        end if
      end do
 
      ! find all polynomials for which exactly one term is quadratic
      ! for these, the Laplacian on the polynomial is 2*metric_inverse(i, i)
      if (powers(2) == 1 .and. powers(0) == dim - 1) then
        do i = 1, dim
          if (polynomials(i, j) == 2) then
            rhs(j) = m_two*metric_inverse(i, i)
          end if
        end do
      end if
      ! find all polynomials for which exactly two terms are linear
      ! for these, the Laplacian on the polynomial is metric_inverse(i, k) + metric_inverse(k, i)
      if (powers(1) == 2 .and. powers(0) == dim - 2) then
        do i = 1, dim
          if (polynomials(i, j) == 1) then
            do k = i+1, dim
              if (polynomials(k, j) == 1) then
                rhs(j) = metric_inverse(i, k) + metric_inverse(k, i)
              end if
            end do
          end if
        end do
      end if
      ! find all polynomials for which exactly one term is linear
      ! for these, the Laplacian on the polynomial is hessian_trace(i)
      ! this term is only non-zero for curvilinear coordinates with a varying metric tensor
      if (powers(1) == 1 .and. powers(0) == dim - 1) then
        do i = 1, dim
          if (polynomials(i, j) == 1) then
            rhs(j) = hessian_trace(i)
          end if
        end do
      end if
    end do
 
    pop_sub(get_rhs_lapl)
  end subroutine get_rhs_lapl
 
  ! ---------------------------------------------------------
  subroutine get_rhs_grad(polynomials, dir, rhs)
    integer, contiguous, intent(in) :: polynomials(:,:)
    integer, intent(in) :: dir
    real(real64), contiguous, intent(out) :: rhs(:)
 
    integer :: j, k, dim
    logical :: this_one
 
    push_sub(get_rhs_grad)
 
    dim = size(polynomials, dim=1)
 
    ! find right-hand side for operator
    rhs(:) = m_zero
    do j = 1, size(polynomials, dim=2)
      this_one = .true.
      do k = 1, dim
        if (k == dir .and. polynomials(k, j) /= 1) this_one = .false.
        if (k /= dir .and. polynomials(k, j) /= 0) this_one = .false.
      end do
      if (this_one) rhs(j) = m_one
    end do
 
    pop_sub(get_rhs_grad)
  end subroutine get_rhs_grad
 
 
  ! ---------------------------------------------------------
  subroutine derivatives_make_discretization(namespace, der, nderiv, ideriv, op, &
    space, name, verbose, order)
    type(namespace_t),           intent(in)    :: namespace
    type(derivatives_t),         intent(in)    :: der
    integer,                     intent(in)    :: nderiv
    integer,                     intent(in)    :: ideriv
    type(nl_operator_t),         intent(inout) :: op(:)
    type(space_t),               intent(in)    :: space
    character(len=80), optional, intent(in)    :: name
    logical, optional,           intent(in)    :: verbose
    integer, optional,           intent(in)    :: order
 
    integer :: p, p_max, i, j, k, pow_max, order_
    real(real64) :: x(der%dim)
    real(real64), allocatable :: mat(:,:), sol(:,:), powers(:,:)
    integer, allocatable :: pol(:,:)
    real(real64), allocatable :: rhs(:,:)
    character(len=80) :: name_
 
    push_sub(derivatives_make_discretization)
 
    order_ = optional_default(order, der%order)
 
    safe_allocate(rhs(1:op(1)%stencil%size, nderiv))
    safe_allocate(pol(1:der%dim, 1:op(1)%stencil%size))
    ! get polynomials
    if (nderiv == 1) then
      if (ideriv <= der%dim) then  ! gradient
        call stencil_pol_grad(der%stencil_type, der%dim, ideriv, order_, pol)
        name_ = index2axis(ideriv) // "-gradient"
      else                      ! Laplacian
        call stencil_pol_lapl(der%stencil_type, op(1)%stencil, der%dim, order_, pol)
        name_ = "Laplacian"
      end if
    else if (nderiv == der%dim) then
      ! in this case, we compute the weights for the gradients at once
      ! because they have the same stencils
      ! this is only used for the cube stencil
      name_ = "gradients"
      call stencil_pol_grad(der%stencil_type, der%dim, 1, order_, pol)
    else
      message(1) = "Error: derivatives_make_discretization can only be called with nderiv = 1 or der%dim."
      call messages_fatal(1)
    end if
 
    safe_allocate(mat(1:op(1)%stencil%size, 1:op(1)%stencil%size))
    safe_allocate(sol(1:op(1)%stencil%size, 1:nderiv))
 
    if (optional_default(verbose, .true.)) then
      name_ = optional_default(name, name_)
      message(1) = 'Info: Generating weights for finite-difference discretization of ' // trim(name_)
      call messages_info(1, namespace=namespace)
    end if
 
    ! use to generate power lookup table
    pow_max = maxval(pol)
    safe_allocate(powers(1:der%dim, 0:pow_max))
    powers(:,:) = m_zero
    powers(:,0) = m_one
 
    p_max = op(1)%np
    if (op(1)%const_w) p_max = 1
 
    do p = 1, p_max
      ! get polynomials and right-hand side
      ! we need the current position for the right-hand side for curvilinear coordinates
      if (nderiv == 1) then
        if (ideriv <= der%dim) then  ! gradient
          call get_rhs_grad(pol, ideriv, rhs(:,1))
          name_ = index2axis(ideriv) // "-gradient"
        else                      ! Laplacian
          call get_rhs_lapl(der%mesh%coord_system, pol, der%mesh%chi(:, p), rhs(:,1), der%stencil_primitive_coordinates)
        end if
      else if (nderiv == der%dim) then
        ! in this case, we compute the weights for the gradients at once
        ! because they have the same stencils
        do i = 1, der%dim
          call get_rhs_grad(pol, i, rhs(:,i))
        end do
      end if
 
      ! first polynomial is just a constant
      mat(1,:) = m_one
      ! i indexes the point in the stencil
      do i = 1, op(1)%stencil%size
        if (ideriv <= der%dim) then
          ! for gradients, we always need to evaluate the polynomials in primitive coordinates
          x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
        else
          ! for the Laplacian, we can evaluate the polynomials in primitive or Cartesian coordinates
          if (der%stencil_primitive_coordinates) then
            x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
          else
            x = der%mesh%x(:, p + op(1)%ri(i, op(1)%rimap(p))) - der%mesh%x(:, p)
          end if
        end if
 
        ! NB: these masses are applied on the cartesian directions. Should add a check for non-orthogonal axes
        x = x*sqrt(der%masses)
 
        ! calculate powers
        powers(:, 1) = x
        do k = 2, pow_max
          powers(:, k) = x*powers(:, k-1)
        end do
 
        ! generate the matrix
        ! j indexes the polynomial being used
        do j = 2, op(1)%stencil%size
          mat(j, i) = powers(1, pol(1, j))
          do k = 2, der%dim
            mat(j, i) = mat(j, i)*powers(k, pol(k, j))
          end do
        end do
      end do ! loop over i = point in stencil
 
      ! linear problem to solve for derivative weights:
      !   mat * sol = rhs
      call lalg_linsyssolve(op(1)%stencil%size, nderiv, mat, rhs, sol)
 
      ! for the cube stencil, all derivatives are calculated at once, so assign
      ! the correct solution to each operator
      do i = 1, nderiv
        op(i)%w(:, p) = sol(:, i)
        ! The central weight of an antisymmetric operator is zero by symmetry
        if (op(i)%symmetry == op_antisymmetric .and. op(i)%stencil%center > 0) then
          op(i)%w(op(i)%stencil%center, p) = m_zero
        endif
      end do
 
    end do ! loop over points p
 
    do i = 1, nderiv
      if (der%remove_zero_weight_points) then
        call nl_operator_remove_zero_weight_points(op(i), space, der%mesh)
      end if
      call nl_operator_output_weights(op(i))
      if (op(i)%const_w) then
        call nl_operator_build_symmetric_weights(op(i))
      end if
    end do
 
    ! In case of constant weights, we store the weights of the Laplacian on the GPU, as this
    ! saves many unecessary transfers
    do i = 1, nderiv
      call nl_operator_allocate_gpu_buffers(op(i))
      call nl_operator_update_gpu_buffers(op(i))
    end do
 
 
    safe_deallocate_a(mat)
    safe_deallocate_a(sol)
    safe_deallocate_a(powers)
    safe_deallocate_a(pol)
    safe_deallocate_a(rhs)
 
    pop_sub(derivatives_make_discretization)
  end subroutine derivatives_make_discretization
 
#ifdef HAVE_MPI
  ! ---------------------------------------------------------
  logical function derivatives_overlap(this) result(overlap)
    type(derivatives_t), intent(in) :: this
 
    push_sub(derivatives_overlap)
 
    overlap = this%comm_method /= blocking
 
    pop_sub(derivatives_overlap)
  end function derivatives_overlap
#endif
 
  ! ---------------------------------------------------------
  subroutine derivatives_get_lapl(this, namespace, op, space, name, order)
    type(derivatives_t), intent(in) :: this
    type(namespace_t), intent(in) :: namespace
    type(nl_operator_t), contiguous, intent(inout) :: op(:)
    class(space_t), intent(in) :: space
    character(len=80), intent(in) :: name
    integer, intent(in) :: order
 
    push_sub(derivatives_get_lapl)
 
    call derivatives_get_stencil_lapl(this, op(1), space, this%mesh%coord_system, name, order)
    call nl_operator_build(space, this%mesh, op(1), this%mesh%np, const_w = .not. this%mesh%use_curvilinear)
    call derivatives_make_discretization(namespace, this, 1, this%dim+1, op(1:1), space, name, order=order)
 
    pop_sub(derivatives_get_lapl)
  end subroutine derivatives_get_lapl
 
  ! ---------------------------------------------------------
  function derivatives_get_inner_boundary_mask(this) result(mask)
    type(derivatives_t), intent(in)  :: this
 
    logical :: mask(1:this%mesh%np)
    integer :: ip, is, index
 
    mask = .false.
    ! Loop through all points in the grid
    do ip = 1, this%mesh%np
      ! For each of them, loop through all points in the stencil
      do is = 1, this%lapl%stencil%size
        ! Get the index of the point obtained as: grid_point + displament_due_to_stencil
        index = nl_operator_get_index(this%lapl, is, ip)
        ! Check whether the displaced point if outsude the grid. Is so, it belongs to the mask
        if (index > this%mesh%np + this%mesh%pv%np_ghost) then
          mask(ip) = .true.
          exit
        end if
      end do
    end do
 
  end function derivatives_get_inner_boundary_mask
 
  ! ---------------------------------------------------------
  real(real64) function derivatives_lapl_get_max_eigenvalue(this)
    type(derivatives_t), intent(in) :: this
 
    integer :: i
 
    push_sub(derivatives_lapl_get_max_eigenvalue)
 
    derivatives_lapl_get_max_eigenvalue = m_zero
    if ((this%stencil_type == der_star .or. this%stencil_type == der_stargeneral) &
      .and. .not. this%mesh%use_curvilinear) then
      ! use Fourier transform of stencil evaluated at the maximum phase
      do i = 1, this%lapl%stencil%size
        derivatives_lapl_get_max_eigenvalue = derivatives_lapl_get_max_eigenvalue + &
          (-1)**maxval(abs(this%lapl%stencil%points(:, i)))*this%lapl%w(i, 1)
      end do
      derivatives_lapl_get_max_eigenvalue = abs(derivatives_lapl_get_max_eigenvalue)
    else
      ! use upper bound from continuum for other stencils
      do i = 1, this%dim
        derivatives_lapl_get_max_eigenvalue = derivatives_lapl_get_max_eigenvalue + &
          m_pi**2/this%mesh%spacing(i)**2
      end do
    end if
 
    pop_sub(derivatives_lapl_get_max_eigenvalue)
  end function derivatives_lapl_get_max_eigenvalue
 
  subroutine derivates_set_coordinates_system(der, coord_system)
    type(derivatives_t),        target, intent(inout) :: der
    class(coordinate_system_t), target, intent(in)    :: coord_system
 
    der%coord_system => coord_system
  end subroutine derivates_set_coordinates_system
 
#include "undef.F90"
#include "real.F90"
#include "derivatives_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "derivatives_inc.F90"
 
end module derivatives_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: