fast3d.f90

Go to the documentation of this file.

! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
!
! The UChicago Argonne, LLC as Operator of Argonne National
! Laboratory holds copyright in the Software. The copyright holder
! reserves all rights except those expressly granted to licensees,
! and U.S. Government license rights.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the disclaimer below.
!
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the disclaimer (as noted below)
! in the documentation and/or other materials provided with the
! distribution.
!
! 3. Neither the name of ANL nor the names of its contributors
! may be used to endorse or promote products derived from this software
! without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
! ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! Additional BSD Notice
! ---------------------
! 1. This notice is required to be provided under our contract with
! the U.S. Department of Energy (DOE). This work was produced at
! Argonne National Laboratory under Contract
! No. DE-AC02-06CH11357 with the DOE.
!
! 2. Neither the United States Government nor UCHICAGO ARGONNE,
! LLC nor any of their employees, makes any warranty,
! express or implied, or assumes any liability or responsibility for the
! accuracy, completeness, or usefulness of any information, apparatus,
! product, or process disclosed, or represents that its use would not
! infringe privately-owned rights.
!
! 3. Also, reference herein to any specific commercial products, process,
! or services by trade name, trademark, manufacturer or otherwise does
! not necessarily constitute or imply its endorsement, recommendation,
! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
! The views and opinions of authors expressed
! herein do not necessarily state or reflect those of the United States
! Government or UCHICAGO ARGONNE, LLC, and shall
! not be used for advertising or product endorsement purposes.
!
module fast3d
  use num_types, only : rp, xp
  use speclib, only : zwgll, zwgl
  use math, only : rzero
  implicit none
  private
 
  public :: fd_weights_full, semhat, setup_intp
 
contains
 
  subroutine fd_weights_full(xi, x, n, m, c)
    integer, intent(in) :: n
    integer, intent(in) :: m
    real(kind=rp), intent(in) :: x(0:n)
    real(kind=rp), intent(out) :: c(0:n,0:m)
    real(kind=rp), intent(in) :: xi
    real(kind=xp) :: cx(0:n,0:m)
    real(kind=xp) :: c1, c2, c3, c4, c5
    integer :: i, j, k, mn
 
    c1 = 1d0
    c4 = x(0) - xi
 
    do k = 0, m
       do j = 0, n
          cx(j,k) = 0.0_xp
       end do
    end do
 
    cx(0,0) = 1d0
 
    do i = 1, n
       mn = min(i,m)
       c2 = 1d0
       c5 = c4
       c4 = x(i) - xi
       do j = 0, i - 1
          c3 = x(i) - x(j)
          c2 = c2 * c3
          do k = mn, 1, -1
             cx(i,k) = c1 * (k * cx(i-1,k-1) - c5 * cx(i-1,k)) / c2
          end do
          cx(i,0) = -c1 * c5 * cx(i-1,0) / c2
          do k = mn, 1, -1
             cx(j,k) = (c4 * cx(j,k) - k * cx(j,k-1)) / c3
          end do
          cx(j,0) = c4 * cx(j,0) / c3
       end do
       c1 = c2
    end do
    c = real(cx, kind=rp)
  end subroutine fd_weights_full
 
 
  subroutine semhat(a, b, c, d, z, dgll, jgll, bgl, zgl, dgl, jgl, n, w)
    integer, intent(in) :: n
    real(kind=rp), intent(inout) :: a(0:n,0:n)
    real(kind=rp), intent(inout) :: b(0:n)
    real(kind=rp), intent(inout) :: c(0:n,0:n)
    real(kind=rp), intent(inout) :: d(0:n,0:n)
    real(kind=rp), intent(inout) :: z(0:n)
    real(kind=rp), intent(inout) :: dgll(0:n,1:n-1), jgll(0:n,1:n-1)
    real(kind=rp), intent(inout) :: bgl(1:n-1)
    real(kind=rp), intent(inout) :: zgl(1:n-1)
    real(kind=rp), intent(inout) :: dgl(1:n-1,0:n)
    real(kind=rp), intent(inout) :: jgl(1:n-1,0:n)
    real(kind=rp), intent(inout) :: w(0:2*n+1)
    integer :: np, nm, n2, i, j, k
    np = n+1
    nm = n-1
    n2 = n-2
    call zwgll(z, b, np)
    do i = 0,n
       call fd_weights_full(z(i), z, n, 1, w)
       do j = 0,n
          d(i,j) = w(j+np) !  Derivative matrix
       end do
    end do
 
    if (n .eq. 1) return ! No interpolation for n = 1
 
    do i = 0,n
       call fd_weights_full(z(i), z(1), n2, 1, w(1))
       do j = 1, nm
          jgll(i,j) = w(j) ! Interpolation matrix
          dgll(i,j) = w(j + nm) ! Derivative matrix
       end do
    end do
    call rzero(a, np*np)
    do j = 0,n
       do i = 0,n
          do k = 0,n
             a(i,j) = a(i,j) + d(k,i)*b(k)*d(k,j)
          end do
          c(i,j) = b(i)*d(i,j)
       end do
    end do
    call zwgl(zgl, bgl, nm)
    do i = 1,n-1
       call fd_weights_full(zgl(i), z, n, 1, w)
       do j = 0,n
          jgl(i,j) = w(j ) !  Interpolation matrix
          dgl(i,j) = w(j+np) !  Derivative    matrix
       end do
    end do
  end subroutine semhat
 
  subroutine setup_intp(jh, jht, z_to, z_from, n_to, n_from, derivative)
    implicit none
    integer, intent(in) :: n_to, n_from, derivative
    real(kind=rp), intent(inout) :: jh(n_to, n_from), jht(n_from, n_to)
    real(kind=rp), intent(in) :: z_to(n_to), z_from(n_from)
    real(kind=rp) :: w(n_from, 0:derivative)
    integer :: i, j
    do i = 1, n_to
       ! This will assign w the weights for interpolating to point
       ! zf(i) values at points zc
       call fd_weights_full(z_to(i), z_from, n_from-1, derivative, w)
 
       ! store each of the weights in the corresponding row/column
       do j = 1, n_from
          jh(i,j) = w(j, derivative)
          jht(j,i) = w(j, derivative)
       end do
    end do
  end subroutine setup_intp
end module fast3d