!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2021 Grzegorz Kowal !! !! 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 3 of the License, 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, see . !! !!***************************************************************************** !! !! module: OPERATORS !! !! This module provides differential operators like gradient, divergence, or !! curl. !! !!***************************************************************************** ! module operators implicit none ! interfaces for procedure pointers ! abstract interface subroutine derivative_1st_iface(d, h, u, v) integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v end subroutine subroutine derivative_2nd_iface(d, h, u, v) integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v end subroutine end interface ! procedure pointers for derivatives ! procedure(derivative_1st_iface), pointer, save :: derivative_1st => null() procedure(derivative_2nd_iface), pointer, save :: derivative_2nd => null() private public :: initialize_operators, finalize_operators public :: divergence, gradient, curl, laplace public :: derivative_1st, derivative_2nd !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_OPERATORS: ! ------------------------------- ! ! Subroutine initializes the module structures, pointers and variables. ! ! Arguments: ! ! order - the order of the interpolation method; ! status - return flag of the procedure execution status; ! !=============================================================================== ! subroutine initialize_operators(order, status) implicit none integer, intent(in) :: order integer, intent(out) :: status !------------------------------------------------------------------------------- ! status = 0 select case(order) case(9) derivative_1st => derivative_1st_9o derivative_2nd => derivative_2nd_9o case(7) derivative_1st => derivative_1st_7o derivative_2nd => derivative_2nd_7o case(5) derivative_1st => derivative_1st_5o derivative_2nd => derivative_2nd_5o case default derivative_1st => derivative_1st_3o derivative_2nd => derivative_2nd_3o end select !------------------------------------------------------------------------------- ! end subroutine initialize_operators ! !=============================================================================== ! ! subroutine FINALIZE_OPERATORS: ! ----------------------------- ! ! Subroutine releases the memory used by module variables and pointers. ! ! Arguments: ! ! status - return flag of the procedure execution status; ! !=============================================================================== ! subroutine finalize_operators(status) implicit none integer, intent(out) :: status !------------------------------------------------------------------------------- ! status = 0 !------------------------------------------------------------------------------- ! end subroutine finalize_operators ! !=============================================================================== ! ! subroutine DIVERGENCE: ! --------------------- ! ! Subroutine calculates the cell centered divergence of the input vector ! field. ! ! div(U) = ∇.([Ux, Uy, Uz]) = ∂x Ux + ∂y Uy + ∂z Uz ! ! Arguments: ! ! dh - the spacial intervals in all direction; ! u - the input vector field; ! v - the output divergence field; ! !=============================================================================== ! subroutine divergence(dh, u, v) implicit none real(kind=8), dimension(3) , intent(in) :: dh real(kind=8), dimension(:,:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:) , intent(out) :: v integer :: dir real(kind=8), dimension(size(u,2), size(u,3), size(u,4)) :: w !------------------------------------------------------------------------------- ! v(:,:,:) = 0.0d+00 ! iterate over directions and update divergence with directional derivatives ! do dir = 1, NDIMS ! calculate derivative along the current direction ! call derivative_1st(dir, dh(dir), u(dir,:,:,:), w(:,:,:)) ! update the divergence array ! v(:,:,:) = v(:,:,:) + w(:,:,:) end do ! directions !------------------------------------------------------------------------------- ! end subroutine divergence ! !=============================================================================== ! ! subroutine GRADIENT: ! ------------------- ! ! Subroutine calculates the cell centered gradient of the input scalar field. ! ! grad(U) = ∇ U = [ ∂x U, ∂y U, ∂z U ] ! ! Arguments: ! ! dh - the spacial intervals in all direction; ! u - the input scalar field; ! v - the output gradient vector field; ! !=============================================================================== ! subroutine gradient(dh, u, v) implicit none real(kind=8), dimension(3) , intent(in) :: dh real(kind=8), dimension(:,:,:) , intent(in) :: u real(kind=8), dimension(:,:,:,:), intent(out) :: v integer :: dir !------------------------------------------------------------------------------- ! v(:,:,:,:) = 0.0d+00 ! iterate over directions and calculate gradient components ! do dir = 1, NDIMS ! calculate derivative along the current direction ! call derivative_1st(dir, dh(dir), u(:,:,:), v(dir,:,:,:)) end do ! directions !------------------------------------------------------------------------------- ! end subroutine gradient ! !=============================================================================== ! ! subroutine CURL: ! --------------- ! ! Subroutine calculates the cell centered curl of the input vector field. ! ! curl(U) = ∇x([Ux, Uy, Uz]) ! = [∂y Uz - ∂z Uy, ∂z Ux - ∂x Uz, ∂x Uy - ∂y Ux] ! ! Arguments: ! ! dh - the spacial intervals in all direction; ! u - the input vector field; ! v - the output divergence field; ! !=============================================================================== ! subroutine curl(dh, u, v) implicit none real(kind=8), dimension(3) , intent(in) :: dh real(kind=8), dimension(:,:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:,:), intent(out) :: v real(kind=8), dimension(size(u,2), size(u,3), size(u,4)) :: w !------------------------------------------------------------------------------- ! ! === calculate Vx component === ! ! contribution from the Y derivative of Uz ! call derivative_1st(2, dh(2), u(3,:,:,:), w(:,:,:)) ! update Vx ! v(1,:,:,:) = w(:,:,:) #if NDIMS == 3 ! contribution from the Z derivative of Uy ! call derivative_1st(3, dh(3), u(2,:,:,:), w(:,:,:)) ! update Vx ! v(1,:,:,:) = v(1,:,:,:) - w(:,:,:) #endif /* NDIMS == 3 */ ! === calculate Vy component === ! #if NDIMS == 3 ! contribution from the Z derivative of Ux ! call derivative_1st(3, dh(3), u(1,:,:,:), w(:,:,:)) ! update Vy ! v(2,:,:,:) = w(:,:,:) ! contribution from the X derivative of Uz ! call derivative_1st(1, dh(1), u(3,:,:,:), w(:,:,:)) ! update Vy ! v(2,:,:,:) = v(2,:,:,:) - w(:,:,:) #else /* NDIMS == 3 */ ! contribution from the X derivative of Az ! call derivative_1st(1, dh(1), u(3,:,:,:), w(:,:,:)) ! update Vy ! v(2,:,:,:) = - w(:,:,:) #endif /* NDIMS == 3 */ ! === calculate Vz component === ! ! contribution from the X derivative of Uy ! call derivative_1st(1, dh(1), u(2,:,:,:), w(:,:,:)) ! update Vz ! v(3,:,:,:) = w(:,:,:) ! contribution from the Y derivative of Ux ! call derivative_1st(2, dh(2), u(1,:,:,:), w(:,:,:)) ! update Vz ! v(3,:,:,:) = v(3,:,:,:) - w(:,:,:) !------------------------------------------------------------------------------- ! end subroutine curl ! !=============================================================================== ! ! subroutine LAPLACE: ! ------------------ ! ! Subroutine calculates the Laplace operator of the input scalar field. ! ! laplace(U) = Δ(U) = ∂²x U + ∂²y U + ∂²z U ! ! Arguments: ! ! dh - the spacial intervals in all direction; ! u - the input scalar field U; ! v - the output field representing the laplacian of u; ! !=============================================================================== ! subroutine laplace(dh, u, v) implicit none real(kind=8), dimension(3) , intent(in) :: dh real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v real(kind=8), dimension(size(u,1), size(u,2), size(u,3)) :: w !------------------------------------------------------------------------------- ! ! calculate the second derivative of U along the X direction ! call derivative_2nd(1, dh(1), u(:,:,:), w(:,:,:)) ! update V ! v(:,:,:) = w(:,:,:) ! calculate the second derivative of U along the X direction ! call derivative_2nd(2, dh(2), u(:,:,:), w(:,:,:)) ! update V ! v(:,:,:) = v(:,:,:) + w(:,:,:) #if NDIMS == 3 ! calculate the second derivative of U along the Z direction ! call derivative_2nd(3, dh(3), u(:,:,:), w(:,:,:)) ! update V ! v(:,:,:) = v(:,:,:) + w(:,:,:) #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine laplace ! !=============================================================================== ! ! subroutine DERIVATIVE_1ST_3O: ! ---------------------------- ! ! Subroutine calculates the first order derivative of the input scalar field ! along the given direction with the 3rd order approximation. ! ! Arguments: ! ! d - the direction of derivative; ! h - the spacial interval; ! u - the input scalar field; ! v - the first derivative of the input field; ! !=============================================================================== ! subroutine derivative_1st_3o(d, h, u, v) implicit none integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v integer :: m0, m1, m2 !------------------------------------------------------------------------------- ! m0 = size(u, d) m1 = m0 - 1 m2 = m0 - 2 select case(d) case(1) v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h v(2:m1,:,:) = 5.0d-01 * (u(3:m0,:,:) - u(1:m2,:,:)) / h v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h case(2) v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h v(:,2:m1,:) = 5.0d-01 * (u(:,3:m0,:) - u(:,1:m2,:)) / h v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h #if NDIMS == 3 case(3) v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h v(:,:,2:m1) = 5.0d-01 * (u(:,:,3:m0) - u(:,:,1:m2)) / h v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine derivative_1st_3o ! !=============================================================================== ! ! subroutine DERIVATIVE_1ST_5O: ! ---------------------------- ! ! Subroutine calculates the first order derivative of the input scalar field ! along the given direction with the 5th order approximation. ! ! Arguments: ! ! d - the direction of derivative; ! h - the spacial interval; ! u - the input scalar field; ! v - the first derivative of the input field; ! !=============================================================================== ! subroutine derivative_1st_5o(d, h, u, v) implicit none integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v integer :: m0, m1, m2, m3, m4 !------------------------------------------------------------------------------- ! m0 = size(u, d) m1 = m0 - 1 m2 = m0 - 2 m3 = m0 - 3 m4 = m0 - 4 select case(d) case(1) v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h v(2 ,:,:) = 5.0d-01 * (u(3 ,:,:) - u(1 ,:,:)) / h v(3:m2,:,:) = (8.0d+00 * (u(4:m1,:,:) - u(2:m3,:,:)) & - (u(5:m0,:,:) - u(1:m4,:,:))) / (1.2d+01 * h) v( m1,:,:) = 5.0d-01 * (u( m0,:,:) - u( m2,:,:)) / h v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h case(2) v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h v(:,2 ,:) = 5.0d-01 * (u(:,3 ,:) - u(:,1 ,:)) / h v(:,3:m2,:) = (8.0d+00 * (u(:,4:m1,:) - u(:,2:m3,:)) & - (u(:,5:m0,:) - u(:,1:m4,:))) / (1.2d+01 * h) v(:, m1,:) = 5.0d-01 * (u(:, m0,:) - u(:, m2,:)) / h v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h #if NDIMS == 3 case(3) v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h v(:,:,2 ) = 5.0d-01 * (u(:,:,3 ) - u(:,:,1 )) / h v(:,:,3:m2) = (8.0d+00 * (u(:,:,4:m1) - u(:,:,2:m3)) & - (u(:,:,5:m0) - u(:,:,1:m4))) / (1.2d+01 * h) v(:,:, m1) = 5.0d-01 * (u(:,:, m0) - u(:,:, m2)) / h v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine derivative_1st_5o ! !=============================================================================== ! ! subroutine DERIVATIVE_1ST_7O: ! ---------------------------- ! ! Subroutine calculates the first order derivative of the input scalar field ! along the given direction with the 7th order approximation. ! ! Arguments: ! ! d - the direction of derivative; ! h - the spacial interval; ! u - the input scalar field; ! v - the first derivative of the input field; ! !=============================================================================== ! subroutine derivative_1st_7o(d, h, u, v) implicit none integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v integer :: m0, m1, m2, m3, m4, m5, m6 !------------------------------------------------------------------------------- ! m0 = size(u, d) m1 = m0 - 1 m2 = m0 - 2 m3 = m0 - 3 m4 = m0 - 4 m5 = m0 - 5 m6 = m0 - 6 select case(d) case(1) v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h v(2 ,:,:) = 5.0d-01 * (u(3 ,:,:) - u(1 ,:,:)) / h v(3 ,:,:) = (8.0d+00 * (u(4 ,:,:) - u(2 ,:,:)) & - (u(5 ,:,:) - u(1 ,:,:))) / (1.2d+01 * h) v(4:m3,:,:) = (4.5d+01 * (u(5:m2,:,:) - u(3:m4,:,:)) & - 9.0d+00 * (u(6:m1,:,:) - u(2:m5,:,:)) & + (u(7:m0,:,:) - u(1:m6,:,:))) / (6.0d+01 * h) v( m2,:,:) = (8.0d+00 * (u( m1,:,:) - u( m3,:,:)) & - (u( m0,:,:) - u( m4,:,:))) / (1.2d+01 * h) v( m1,:,:) = 5.0d-01 * (u( m0,:,:) - u( m2,:,:)) / h v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h case(2) v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h v(:,2 ,:) = 5.0d-01 * (u(:,3 ,:) - u(:,1 ,:)) / h v(:,3 ,:) = (8.0d+00 * (u(:,4 ,:) - u(:,2 ,:)) & - (u(:,5 ,:) - u(:,1 ,:))) / (1.2d+01 * h) v(:,4:m3,:) = (4.5d+01 * (u(:,5:m2,:) - u(:,3:m4,:)) & - 9.0d+00 * (u(:,6:m1,:) - u(:,2:m5,:)) & + (u(:,7:m0,:) - u(:,1:m6,:))) / (6.0d+01 * h) v(:, m2,:) = (8.0d+00 * (u(:, m1,:) - u(:, m3,:)) & - (u(:, m0,:) - u(:, m4,:))) / (1.2d+01 * h) v(:, m1,:) = 5.0d-01 * (u(:, m0,:) - u(:, m2,:)) / h v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h #if NDIMS == 3 case(3) v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h v(:,:,2 ) = 5.0d-01 * (u(:,:,3 ) - u(:,:,1 )) / h v(:,:,3 ) = (8.0d+00 * (u(:,:,4 ) - u(:,:,2 )) & - (u(:,:,5 ) - u(:,:,1 ))) / (1.2d+01 * h) v(:,:,4:m3) = (4.5d+01 * (u(:,:,5:m2) - u(:,:,3:m4)) & - 9.0d+00 * (u(:,:,6:m1) - u(:,:,2:m5)) & + (u(:,:,7:m0) - u(:,:,1:m6))) / (6.0d+01 * h) v(:,:, m2) = (8.0d+00 * (u(:,:, m1) - u(:,:, m3)) & - (u(:,:, m0) - u(:,:, m4))) / (1.2d+01 * h) v(:,:, m1) = 5.0d-01 * (u(:,:, m0) - u(:,:, m2)) / h v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine derivative_1st_7o ! !=============================================================================== ! ! subroutine DERIVATIVE_1ST_9O: ! ---------------------------- ! ! Subroutine calculates the first order derivative of the input scalar field ! along the given direction with the 9th order approximation. ! ! Arguments: ! ! d - the direction of derivative; ! h - the spacial interval; ! u - the input scalar field; ! v - the first derivative of the input field; ! !=============================================================================== ! subroutine derivative_1st_9o(d, h, u, v) implicit none integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v integer :: m0, m1, m2, m3, m4, m5, m6, m7, m8 !------------------------------------------------------------------------------- ! m0 = size(u, d) m1 = m0 - 1 m2 = m0 - 2 m3 = m0 - 3 m4 = m0 - 4 m5 = m0 - 5 m6 = m0 - 6 m7 = m0 - 7 m8 = m0 - 8 select case(d) case(1) v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h v(2 ,:,:) = 5.00d-01 * (u(3 ,:,:) - u(1 ,:,:)) / h v(3 ,:,:) = (8.00d+00 * (u(4 ,:,:) - u(2 ,:,:)) & - (u(5 ,:,:) - u(1 ,:,:))) / (1.2d+01 * h) v(4 ,:,:) = (4.50d+01 * (u(5 ,:,:) - u(3 ,:,:)) & - 9.00d+00 * (u(6 ,:,:) - u(2 ,:,:)) & + (u(7 ,:,:) - u(1 ,:,:))) / (6.0d+01 * h) v(5:m4,:,:) = (6.72d+02 * (u(6:m3,:,:) - u(4:m5,:,:)) & - 1.68d+02 * (u(7:m2,:,:) - u(3:m6,:,:)) & + 3.20d+01 * (u(8:m1,:,:) - u(2:m7,:,:)) & - 3.00d+00 * (u(9:m0,:,:) - u(1:m8,:,:))) / (8.4d+02 * h) v( m3,:,:) = (4.50d+01 * (u( m2,:,:) - u( m4,:,:)) & - 9.00d+00 * (u( m1,:,:) - u( m5,:,:)) & + (u( m0,:,:) - u( m6,:,:))) / (6.0d+01 * h) v( m2,:,:) = (8.00d+00 * (u( m1,:,:) - u( m3,:,:)) & - (u( m0,:,:) - u( m4,:,:))) / (1.2d+01 * h) v( m1,:,:) = 5.00d-01 * (u( m0,:,:) - u( m2,:,:)) / h v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h case(2) v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h v(:,2 ,:) = 5.00d-01 * (u(:,3 ,:) - u(:,1 ,:)) / h v(:,3 ,:) = (8.00d+00 * (u(:,4 ,:) - u(:,2 ,:)) & - (u(:,5 ,:) - u(:,1 ,:))) / (1.2d+01 * h) v(:,4 ,:) = (4.50d+01 * (u(:,5 ,:) - u(:,3 ,:)) & - 9.00d+00 * (u(:,6 ,:) - u(:,2 ,:)) & + (u(:,7 ,:) - u(:,1 ,:))) / (6.0d+01 * h) v(:,5:m4,:) = (6.72d+02 * (u(:,6:m3,:) - u(:,4:m5,:)) & - 1.68d+02 * (u(:,7:m2,:) - u(:,3:m6,:)) & + 3.20d+01 * (u(:,8:m1,:) - u(:,2:m7,:)) & - 3.00d+00 * (u(:,9:m0,:) - u(:,1:m8,:))) / (8.4d+02 * h) v(:, m3,:) = (4.50d+01 * (u(:, m2,:) - u(:, m4,:)) & - 9.00d+00 * (u(:, m1,:) - u(:, m5,:)) & + (u(:, m0,:) - u(:, m6,:))) / (6.0d+01 * h) v(:, m2,:) = (8.00d+00 * (u(:, m1,:) - u(:, m3,:)) & - (u(:, m0,:) - u(:, m4,:))) / (1.2d+01 * h) v(:, m1,:) = 5.00d-01 * (u(:, m0,:) - u(:, m2,:)) / h v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h #if NDIMS == 3 case(3) v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h v(:,:,2 ) = 5.00d-01 * (u(:,:,3 ) - u(:,:,1 )) / h v(:,:,3 ) = (8.00d+00 * (u(:,:,4 ) - u(:,:,2 )) & - (u(:,:,5 ) - u(:,:,1 ))) / (1.2d+01 * h) v(:,:,4 ) = (4.50d+01 * (u(:,:,5 ) - u(:,:,3 )) & - 9.00d+00 * (u(:,:,6 ) - u(:,:,2 )) & + (u(:,:,7 ) - u(:,:,1 ))) / (6.0d+01 * h) v(:,:,5:m4) = (6.72d+02 * (u(:,:,6:m3) - u(:,:,4:m5)) & - 1.68d+02 * (u(:,:,7:m2) - u(:,:,3:m6)) & + 3.20d+01 * (u(:,:,8:m1) - u(:,:,2:m7)) & - 3.00d+00 * (u(:,:,9:m0) - u(:,:,1:m8))) / (8.4d+02 * h) v(:,:, m3) = (4.50d+01 * (u(:,:, m2) - u(:,:, m4)) & - 9.00d+00 * (u(:,:, m1) - u(:,:, m5)) & + (u(:,:, m0) - u(:,:, m6))) / (6.0d+01 * h) v(:,:, m2) = (8.00d+00 * (u(:,:, m1) - u(:,:, m3)) & - (u(:,:, m0) - u(:,:, m4))) / (1.2d+01 * h) v(:,:, m1) = 5.00d-01 * (u(:,:, m0) - u(:,:, m2)) / h v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine derivative_1st_9o ! !=============================================================================== ! ! subroutine DERIVATIVE_2ND_3O: ! ---------------------------- ! ! Subroutine calculates the second order derivative of the input scalar field ! along the given direction with the 3rd order approximation. ! ! Arguments: ! ! d - the direction of derivative; ! h - the spacial interval; ! u - the input scalar field; ! v - the output scalar field representing the second derivative of u; ! !=============================================================================== ! subroutine derivative_2nd_3o(d, h, u, v) implicit none integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v integer :: m0, m1, m2 real(kind=8) :: h2 !------------------------------------------------------------------------------- ! m0 = size(u, d) m1 = m0 - 1 m2 = m0 - 2 h2 = h * h select case(d) case(1) v(1 ,:,:) = 0.0d+00 v(2:m1,:,:) = ((u(3:m0,:,:) + u(1:m2,:,:)) - 2.0d+00 * u(2:m1,:,:)) / h2 v( m0,:,:) = 0.0d+00 case(2) v(:,1 ,:) = 0.0d+00 v(:,2:m1,:) = ((u(:,3:m0,:) + u(:,1:m2,:)) - 2.0d+00 * u(:,2:m1,:)) / h2 v(:, m0,:) = 0.0d+00 #if NDIMS == 3 case(3) v(:,:,1 ) = 0.0d+00 v(:,:,2:m1) = ((u(:,:,3:m0) + u(:,:,1:m2)) - 2.0d+00 * u(:,:,2:m1)) / h2 v(:,:, m0) = 0.0d+00 #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine derivative_2nd_3o ! !=============================================================================== ! ! subroutine DERIVATIVE_2ND_5O: ! ---------------------------- ! ! Subroutine calculates the second order derivative of the input scalar field ! along the given direction with the 5th order approximation. ! ! Arguments: ! ! d - the direction of derivative; ! h - the spacial interval; ! u - the input scalar field; ! v - the output scalar field representing the second derivative of u; ! !=============================================================================== ! subroutine derivative_2nd_5o(d, h, u, v) implicit none integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v integer :: m0, m1, m2, m3, m4 real(kind=8) :: h2 !------------------------------------------------------------------------------- ! m0 = size(u, d) m1 = m0 - 1 m2 = m0 - 2 m3 = m0 - 3 m4 = m0 - 4 h2 = h * h select case(d) case(1) v(1 ,:,:) = 0.0d+00 v(2 ,:,:) = ((u(3 ,:,:) + u(1 ,:,:)) - 2.0d+00 * u(2 ,:,:)) / h2 v(3:m2,:,:) = (1.6d+01 * (u(4:m1,:,:) + u(2:m3,:,:)) & - (u(5:m0,:,:) + u(1:m4,:,:)) & - 3.0d+01 * u(3:m2,:,:)) / (1.2d+01 * h2) v( m1,:,:) = ((u( m0,:,:) + u( m2,:,:)) - 2.0d+00 * u( m1,:,:)) / h2 v( m0,:,:) = 0.0d+00 case(2) v(:,1 ,:) = 0.0d+00 v(:,2 ,:) = ((u(:,3 ,:) + u(:,1 ,:)) - 2.0d+00 * u(:,2 ,:)) / h2 v(:,3:m2,:) = (1.6d+01 * (u(:,4:m1,:) + u(:,2:m3,:)) & - (u(:,5:m0,:) + u(:,1:m4,:)) & - 3.0d+01 * u(:,3:m2,:)) / (1.2d+01 * h2) v(:, m1,:) = ((u(:, m0,:) + u(:, m2,:)) - 2.0d+00 * u(:, m1,:)) / h2 v(:, m0,:) = 0.0d+00 #if NDIMS == 3 case(3) v(:,:,1 ) = 0.0d+00 v(:,:,2 ) = ((u(:,:,3 ) + u(:,:,1 )) - 2.0d+00 * u(:,:,2 )) / h2 v(:,:,3:m2) = (1.6d+01 * (u(:,:,4:m1) + u(:,:,2:m3)) & - (u(:,:,5:m0) + u(:,:,1:m4)) & - 3.0d+01 * u(:,:,3:m2)) / (1.2d+01 * h2) v(:,:, m1) = ((u(:,:, m0) + u(:,:, m2)) - 2.0d+00 * u(:,:, m1)) / h2 v(:,:, m0) = 0.0d+00 #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine derivative_2nd_5o ! !=============================================================================== ! ! subroutine DERIVATIVE_2ND_7O: ! ---------------------------- ! ! Subroutine calculates the second order derivative of the input scalar field ! along the given direction with the 7th order approximation. ! ! Arguments: ! ! d - the direction of derivative; ! h - the spacial interval; ! u - the input scalar field; ! v - the output scalar field representing the second derivative of u; ! !=============================================================================== ! subroutine derivative_2nd_7o(d, h, u, v) implicit none integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v integer :: m0, m1, m2, m3, m4, m5, m6 real(kind=8) :: h2 !------------------------------------------------------------------------------- ! m0 = size(u, d) m1 = m0 - 1 m2 = m0 - 2 m3 = m0 - 3 m4 = m0 - 4 m5 = m0 - 5 m6 = m0 - 6 h2 = h * h select case(d) case(1) v(1 ,:,:) = 0.0d+00 v(2 ,:,:) = ((u(3 ,:,:) + u(1 ,:,:)) - 2.0d+00 * u(2 ,:,:)) / h2 v(3 ,:,:) = (1.60d+01 * (u(4 ,:,:) + u(2 ,:,:)) & - (u(5 ,:,:) + u(1 ,:,:)) & - 3.00d+01 * u(3 ,:,:)) / (1.2d+01 * h2) v(4:m3,:,:) = (1.35d+02 * (u(5:m2,:,:) + u(3:m4,:,:)) & - 1.35d+01 * (u(6:m1,:,:) + u(2:m5,:,:)) & + (u(7:m0,:,:) + u(1:m6,:,:)) & - 2.45d+02 * u(4:m3,:,:)) / (9.0d+01 * h2) v( m2,:,:) = (1.60d+01 * (u( m1,:,:) + u( m3,:,:)) & - (u( m0,:,:) + u( m4,:,:)) & - 3.00d+01 * u( m2,:,:)) / (1.2d+01 * h2) v( m1,:,:) = ((u( m0,:,:) + u( m2,:,:)) - 2.0d+00 * u( m1,:,:)) / h2 v( m0,:,:) = 0.0d+00 case(2) v(:,1 ,:) = 0.0d+00 v(:,2 ,:) = ((u(:,3 ,:) + u(:,1 ,:)) - 2.0d+00 * u(:,2 ,:)) / h2 v(:,3 ,:) = (1.60d+01 * (u(:,4 ,:) + u(:,2 ,:)) & - (u(:,5 ,:) + u(:,1 ,:)) & - 3.00d+01 * u(:,3 ,:)) / (1.2d+01 * h2) v(:,4:m3,:) = (1.35d+02 * (u(:,5:m2,:) + u(:,3:m4,:)) & - 1.35d+01 * (u(:,6:m1,:) + u(:,2:m5,:)) & + (u(:,7:m0,:) + u(:,1:m6,:)) & - 2.45d+02 * u(:,4:m3,:)) / (9.0d+01 * h2) v(:, m2,:) = (1.60d+01 * (u(:, m1,:) + u(:, m3,:)) & - (u(:, m0,:) + u(:, m4,:)) & - 3.00d+01 * u(:, m2,:)) / (1.2d+01 * h2) v(:, m1,:) = ((u(:, m0,:) + u(:, m2,:)) - 2.0d+00 * u(:, m1,:)) / h2 v(:, m0,:) = 0.0d+00 #if NDIMS == 3 case(3) v(:,:,1 ) = 0.0d+00 v(:,:,2 ) = ((u(:,:,3 ) + u(:,:,1 )) - 2.0d+00 * u(:,:,2 )) / h2 v(:,:,3 ) = (1.60d+01 * (u(:,:,4 ) + u(:,:,2 )) & - (u(:,:,5 ) + u(:,:,1 )) & - 3.00d+01 * u(:,:,3 )) / (1.2d+01 * h2) v(:,:,4:m3) = (1.35d+02 * (u(:,:,5:m2) + u(:,:,3:m4)) & - 1.35d+01 * (u(:,:,6:m1) + u(:,:,2:m5)) & + (u(:,:,7:m0) + u(:,:,1:m6)) & - 2.45d+02 * u(:,:,4:m3)) / (9.0d+01 * h2) v(:,:, m2) = (1.60d+01 * (u(:,:, m1) + u(:,:, m3)) & - (u(:,:, m0) + u(:,:, m4)) & - 3.00d+01 * u(:,:, m2)) / (1.2d+01 * h2) v(:,:, m1) = ((u(:,:, m0) + u(:,:, m2)) - 2.0d+00 * u(:,:, m1)) / h2 v(:,:, m0) = 0.0d+00 #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine derivative_2nd_7o ! !=============================================================================== ! ! subroutine DERIVATIVE_2ND_9O: ! ---------------------------- ! ! Subroutine calculates the second order derivative of the input scalar field ! along the given direction with the 9th order approximation. ! ! Arguments: ! ! d - the direction of derivative; ! h - the spacial interval; ! u - the input scalar field; ! v - the output scalar field representing the second derivative of u; ! !=============================================================================== ! subroutine derivative_2nd_9o(d, h, u, v) implicit none integer , intent(in) :: d real(kind=8) , intent(in) :: h real(kind=8), dimension(:,:,:), intent(in) :: u real(kind=8), dimension(:,:,:), intent(out) :: v integer :: m0, m1, m2, m3, m4, m5, m6, m7, m8 real(kind=8) :: h2 !------------------------------------------------------------------------------- ! m0 = size(u, d) m1 = m0 - 1 m2 = m0 - 2 m3 = m0 - 3 m4 = m0 - 4 m5 = m0 - 5 m6 = m0 - 6 m7 = m0 - 7 m8 = m0 - 8 h2 = h * h select case(d) case(1) v(1 ,:,:) = 0.0d+00 v(2 ,:,:) = ((u(3 ,:,:) + u(1 ,:,:)) - 2.0d+00 * u(2 ,:,:)) / h2 v(3 ,:,:) = (1.600d+01 * (u(4 ,:,:) + u(2 ,:,:)) & - (u(5 ,:,:) + u(1 ,:,:)) & - 3.000d+01 * u(3 ,:,:)) / (1.2d+01 * h2) v(4 ,:,:) = (1.350d+02 * (u(5 ,:,:) + u(3 ,:,:)) & - 1.350d+01 * (u(6 ,:,:) + u(2 ,:,:)) & + (u(7 ,:,:) + u(1 ,:,:)) & - 2.450d+02 * u(4 ,:,:)) / (9.0d+01 * h2) v(5:m4,:,:) = (8.064d+03 * (u(6:m3,:,:) + u(4:m5,:,:)) & - 1.008d+03 * (u(7:m2,:,:) + u(3:m6,:,:)) & + 1.280d+02 * (u(8:m1,:,:) + u(2:m7,:,:)) & - 9.000d+00 * (u(9:m0,:,:) + u(1:m8,:,:)) & - 1.435d+04 * u(5:m4,:,:)) / (5.04d+03 * h2) v( m3,:,:) = (1.350d+02 * (u( m2,:,:) + u( m4,:,:)) & - 1.350d+01 * (u( m1,:,:) + u( m5,:,:)) & + (u( m0,:,:) + u( m6,:,:)) & - 2.450d+02 * u( m3,:,:)) / (9.0d+01 * h2) v( m2,:,:) = (1.600d+01 * (u( m1,:,:) + u( m3,:,:)) & - (u( m0,:,:) + u( m4,:,:)) & - 3.000d+01 * u( m2,:,:)) / (1.2d+01 * h2) v( m1,:,:) = ((u( m0,:,:) + u( m2,:,:)) - 2.0d+00 * u( m1,:,:)) / h2 v( m0,:,:) = 0.0d+00 case(2) v(:,1 ,:) = 0.0d+00 v(:,2 ,:) = ((u(:,3 ,:) + u(:,1 ,:)) - 2.0d+00 * u(:,2 ,:)) / h2 v(:,3 ,:) = (1.600d+01 * (u(:,4 ,:) + u(:,2 ,:)) & - (u(:,5 ,:) + u(:,1 ,:)) & - 3.000d+01 * u(:,3 ,:)) / (1.2d+01 * h2) v(:,4 ,:) = (1.350d+02 * (u(:,5 ,:) + u(:,3 ,:)) & - 1.350d+01 * (u(:,6 ,:) + u(:,2 ,:)) & + (u(:,7 ,:) + u(:,1 ,:)) & - 2.450d+02 * u(:,4 ,:)) / (9.0d+01 * h2) v(:,5:m4,:) = (8.064d+03 * (u(:,6:m3,:) + u(:,4:m5,:)) & - 1.008d+03 * (u(:,7:m2,:) + u(:,3:m6,:)) & + 1.280d+02 * (u(:,8:m1,:) + u(:,2:m7,:)) & - 9.000d+00 * (u(:,9:m0,:) + u(:,1:m8,:)) & - 1.435d+04 * u(:,5:m4,:)) / (5.04d+03 * h2) v(:, m3,:) = (1.350d+02 * (u(:, m2,:) + u(:, m4,:)) & - 1.350d+01 * (u(:, m1,:) + u(:, m5,:)) & + (u(:, m0,:) + u(:, m6,:)) & - 2.450d+02 * u(:, m3,:)) / (9.0d+01 * h2) v(:, m2,:) = (1.600d+01 * (u(:, m1,:) + u(:, m3,:)) & - (u(:, m0,:) + u(:, m4,:)) & - 3.000d+01 * u(:, m2,:)) / (1.2d+01 * h2) v(:, m1,:) = ((u(:, m0,:) + u(:, m2,:)) - 2.0d+00 * u(:, m1,:)) / h2 v(:, m0,:) = 0.0d+00 #if NDIMS == 3 case(3) v(:,:,1 ) = 0.0d+00 v(:,:,2 ) = ((u(:,:,3 ) + u(:,:,1 )) - 2.0d+00 * u(:,:,2 )) / h2 v(:,:,3 ) = (1.600d+01 * (u(:,:,4 ) + u(:,:,2 )) & - (u(:,:,5 ) + u(:,:,1 )) & - 3.000d+01 * u(:,:,3 )) / (1.2d+01 * h2) v(:,:,4 ) = (1.350d+02 * (u(:,:,5 ) + u(:,:,3 )) & - 1.350d+01 * (u(:,:,6 ) + u(:,:,2 )) & + (u(:,:,7 ) + u(:,:,1 )) & - 2.450d+02 * u(:,:,4 )) / (9.0d+01 * h2) v(:,:,5:m4) = (8.064d+03 * (u(:,:,6:m3) + u(:,:,4:m5)) & - 1.008d+03 * (u(:,:,7:m2) + u(:,:,3:m6)) & + 1.280d+02 * (u(:,:,8:m1) + u(:,:,2:m7)) & - 9.000d+00 * (u(:,:,9:m0) + u(:,:,1:m8)) & - 1.435d+04 * u(:,:,5:m4)) / (5.04d+03 * h2) v(:,:, m3) = (1.350d+02 * (u(:,:, m2) + u(:,:, m4)) & - 1.350d+01 * (u(:,:, m1) + u(:,:, m5)) & + (u(:,:, m0) + u(:,:, m6)) & - 2.450d+02 * u(:,:, m3)) / (9.0d+01 * h2) v(:,:, m2) = (1.600d+01 * (u(:,:, m1) + u(:,:, m3)) & - (u(:,:, m0) + u(:,:, m4)) & - 3.000d+01 * u(:,:, m2)) / (1.2d+01 * h2) v(:,:, m1) = ((u(:,:, m0) + u(:,:, m2)) - 2.0d+00 * u(:,:, m1)) / h2 v(:,:, m0) = 0.0d+00 #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine derivative_2nd_9o !=============================================================================== ! end module operators