derivatives_test.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_test_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use boundaries_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use lalg_adv_oct_m
  use loct_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use utils_oct_m
 
  implicit none
 
  private
  public ::                             &
    dderivatives_test,                  &
    zderivatives_test,                  &
    derivatives_advanced_benchmark
 
contains
 
  subroutine derivatives_advanced_benchmark(this, namespace)
    type(derivatives_t), intent(in) :: this
    type(namespace_t),   intent(in) :: namespace
 
    real(real64), allocatable :: ff(:), grad(:,:), resgrad(:, :), lapl(:), reslapl(:)
    real(real64) :: kvec(3)
 
    safe_allocate(ff(1:this%mesh%np_part))
    safe_allocate(grad(1:this%mesh%np, 1:this%mesh%box%dim))
    safe_allocate(lapl(1:this%mesh%np))
    safe_allocate(reslapl(1:this%mesh%np))
    safe_allocate(resgrad(1:this%mesh%np, 1:this%mesh%box%dim))
 
    ! Testing on some planewaves
    kvec = [0.5_real64, 0.0_real64, 0.0_real64]         ! low frequency
    call testfunc_plane_wave(this%mesh, ff, grad, lapl, kvec)
    call test_function("Testing Planewave k=(0.5,0.0,0.0)")
 
    kvec = [1.0_real64, 1.0_real64, 0.0_real64]         ! diagonal
    call testfunc_plane_wave(this%mesh, ff, grad, lapl, kvec)
    call test_function("Testing Planewave k=(1.0,1.0,0.0)")
 
    kvec = [m_pi/this%mesh%spacing(1), 0.0_real64, 0.0_real64]            ! Nyquist
    call testfunc_plane_wave(this%mesh, ff, grad, lapl, kvec)
    call test_function("Testing Planewave k=(\pi/dr,0.0,0.0)")
 
    kvec = [m_pi/this%mesh%spacing(1), m_pi/this%mesh%spacing(2), 0.0_real64]         ! diagonal Nyquist
    call testfunc_plane_wave(this%mesh, ff, grad, lapl, kvec)
    call test_function("Testing Planewave k=(\pi/dr,\pi/dr,0.0)")
 
    ! Testing on some polynomials
    call testfunc_polynomial(this%mesh, ff, grad, lapl, p=2)
    call test_function("Testing degree-2 polynomial")
 
    call testfunc_polynomial(this%mesh, ff, grad, lapl, p=3)
    call test_function("Testing degree-3 polynomial")
 
    call testfunc_polynomial(this%mesh, ff, grad, lapl, p=4)
    call test_function("Testing degree-4 polynomial")
 
    call testfunc_polynomial(this%mesh, ff, grad, lapl, p=6)
    call test_function("Testing degree-6 polynomial")
 
    call testfunc_polynomial(this%mesh, ff, grad, lapl, p=8)
    call test_function("Testing degree-8 polynomial")
 
 
    safe_deallocate_a(ff)
    safe_deallocate_a(grad)
    safe_deallocate_a(lapl)
    safe_deallocate_a(resgrad)
    safe_deallocate_a(reslapl)
 
  contains
    subroutine test_function(label)
      character(len=*), intent(in) :: label
 
      message(1) = trim(label)
 
      call dderivatives_grad(this, ff, resgrad, set_bc=.false.)
      resgrad = resgrad - grad
      write(message(2), '(a, es17.10)') 'Gradient err = ', dmf_nrm2(this%mesh, this%mesh%box%dim, resgrad)
 
      call dderivatives_lapl(this, ff, reslapl, set_bc=.false.)
      reslapl = reslapl - lapl
      write(message(3), '(a, es17.10)') 'Laplacian err = ', dmf_nrm2(this%mesh, reslapl)
 
      call messages_info(3, namespace=namespace)
    end subroutine test_function
  end subroutine derivatives_advanced_benchmark
 
  subroutine testfunc_plane_wave(mesh, ff, grad, lapl, kvec)
    type(mesh_t), intent(in) :: mesh
    real(real64), intent(out) :: ff(:), grad(:, :), lapl(:)
    real(real64), intent(in) :: kvec(3)
 
    integer :: ip, dim, idir
    real(real64) :: phase, kdotx, kk
 
    dim = mesh%box%dim
 
    kk = dot_product(kvec(1:dim), kvec(1:dim))
    do ip = 1, mesh%np
      kdotx = dot_product(kvec(1:dim), mesh%x(1:dim, ip))
      phase = cos(kdotx)
 
      ff(ip) = phase
 
      ! exact gradient
      do idir = 1, mesh%box%dim
        grad(ip, idir) = -kvec(idir)*sin(kdotx)
      end do
 
      ! exact Laplacian
      lapl(ip) = -kk * phase
    end do
 
    do ip = mesh%np+1, mesh%np_part
      kdotx = dot_product(kvec(1:dim), mesh%x(1:dim, ip))
      phase = cos(kdotx)
 
      ff(ip) = phase
    end do
  end subroutine testfunc_plane_wave
 
  subroutine testfunc_polynomial(mesh, ff, grad, lapl, p)
    type(mesh_t), intent(in) :: mesh
    real(real64), intent(out) :: ff(:), grad(:, :), lapl(:)
    integer, intent(in) :: p
    integer :: ip, d, dim
    real(real64) :: x(3), r2
 
    dim = mesh%box%dim
 
    do ip = 1, mesh%np
      x = 0.0_real64
      x(1:dim) = mesh%x(1:dim, ip)
 
      ! f = (x^2 + y^2 + z^2)^(p/2)
      r2 = dot_product(x(1:dim), x(1:dim))
      ff(ip) = r2**(0.5_real64*p)
 
      ! gradient: grad(f) = p * r^(p-2) * x
      if (p >= 2) then
        do d = 1, dim
          grad(ip,d) = p * r2**(0.5_real64*p - 1.0_real64) * x(d)
        end do
      else
        grad(ip,:) = 0.0_real64
      end if
 
      ! Laplacian: ∇² r^n = n (n + d - 2) r^(n-2)
      if (p >= 2) then
        lapl(ip) = p * (p + dim - 2) * r2**(0.5_real64*p - 1.0_real64)
      else
        lapl(ip) = 0.0_real64
      end if
    end do
 
    do ip = mesh%np+1, mesh%np_part
      x = 0.0_real64
      x(1:dim) = mesh%x(1:dim, ip)
 
      ! f = (x^2 + y^2 + z^2)^(p/2)
      r2 = dot_product(x(1:dim), x(1:dim))
      ff(ip) = r2**(0.5_real64*p)
    end do
  end subroutine testfunc_polynomial
 
 
#include "undef.F90"
#include "real.F90"
#include "derivatives_test_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "derivatives_test_inc.F90"
 
end module derivatives_test_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: