stencil_variational.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 stencil_variational_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use messages_oct_m
  use profiling_oct_m
  use nl_operator_oct_m
 
  implicit none
 
  private
  public ::                         &
    stencil_variational_coeff_lapl
 
contains
 
  ! ---------------------------------------------------------
  subroutine stencil_variational_coeff_lapl(dim, order, h, lapl, alpha)
    integer,             intent(in)    :: dim
    integer,             intent(in)    :: order
    real(real64),        intent(in)    :: h(:)
    type(nl_operator_t), intent(inout) :: lapl
    real(real64), optional,     intent(in)    :: alpha
 
    integer :: i, j, k
    real(real64) :: alpha_, kmax
    real(real64), allocatable :: fp(:)
 
    push_sub(stencil_variational_coeff_lapl)
 
    alpha_ = m_one
    if (present(alpha)) alpha_ = alpha
    kmax = m_pi**2 * alpha_
 
    safe_allocate(fp(1:order+1))
    select case (order)
    case (1)
      fp(1) = -kmax/m_two
      fp(2) =  kmax/m_four
    case (2)
      fp(1) = -m_half-m_three*kmax/8.0_real64
      fp(2) =  kmax/m_four
      fp(3) =  m_one/m_four - kmax/16.0_real64
    case (3)
      fp(1) = -m_five/6.0_real64 - m_five*kmax/16.0_real64
      fp(2) =  m_one/12.0_real64 + 15.0_real64*kmax/64.0_real64
      fp(3) =  m_five/12.0_real64 - m_three*kmax/32.0_real64
      fp(4) = -m_one/12.0_real64 + kmax/64.0_real64
    case (4)
      fp(1) = -77.0_real64/72.0_real64 - 35.0_real64*kmax/128.0_real64
      fp(2) =  8.0_real64/45.0_real64 + 7.0_real64*kmax/32.0_real64
      fp(3) =  23.0_real64/45.0_real64 - 7.0_real64*kmax/64.0_real64
      fp(4) = -8.0_real64/45.0_real64 + kmax/32.0_real64
      fp(5) =  17.0_real64/720.0_real64 - kmax/256.0_real64
    case (5)
      fp(1) = -449.0_real64/360.0_real64 - 63.0_real64*kmax/256.0_real64
      fp(2) =  m_four/15.0_real64 + 105.0_real64*kmax/512.0_real64
      fp(3) =  59.0_real64/105.0_real64 - 15.0_real64*kmax/128.0_real64
      fp(4) = -82.0_real64/315.0_real64 + 45.0_real64*kmax/1024.0_real64
      fp(5) =  311.0_real64/5040.0_real64 - m_five*kmax/512.0_real64
      fp(6) = -m_two/315.0_real64 + kmax/1024.0_real64
    case (6)
      fp(1) = -2497.0_real64/1800.0_real64 - 231.0_real64*kmax/1024.0_real64
      fp(2) =  26.0_real64/75.0_real64 + 99.0_real64*kmax/512.0_real64
      fp(3) =  493.0_real64/840.0_real64 - 495.0_real64*kmax/4096.0_real64
      fp(4) = -103.0_real64/315.0_real64 + 55.0_real64*kmax/1024.0_real64
      fp(5) =  2647.0_real64/25200.0_real64 - 33.0_real64*kmax/2048.0_real64
      fp(6) = -31.0_real64/1575.0_real64 + m_three*kmax/1024.0_real64
      fp(7) =  m_one/600.0_real64 - kmax/4096.0_real64
    end select
 
    lapl%w(1,:) = fp(1)*sum(1/h(1:dim)**2)
 
    k = 1
    do i = 1, dim
      do j = -order, -1
        k = k + 1
        lapl%w(k,:) = fp(-j+1) / h(i)**2
      end do
 
      do j = 1, order
        k = k + 1
        lapl%w(k,:) = fp( j+1) / h(i)**2
      end do
    end do
 
    safe_deallocate_a(fp)
 
    pop_sub(stencil_variational_coeff_lapl)
  end subroutine stencil_variational_coeff_lapl
 
 
end module stencil_variational_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: