From 0c913f2c81ce7f1b10c89e38299badf23747763e Mon Sep 17 00:00:00 2001 From: diener Date: Mon, 10 Apr 2006 15:10:43 +0000 Subject: First attempt at upwinding SBP operators. Don't use yet. They will change soon to be more flexible. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@62 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- src/Derivatives_2_1.F90 | 183 ++++++++++ src/Derivatives_4_2.F90 | 292 ++++++++++++++++ src/Derivatives_6_3_min_err_coeff.F90 | 426 +++++++++++++++++++++++ src/Derivatives_8_4_min_err_coeff.F90 | 619 ++++++++++++++++++++++++++++++++++ src/call_up_derivs.c | 206 +++++++++++ src/make.code.defn | 1 + 6 files changed, 1727 insertions(+) create mode 100644 src/call_up_derivs.c diff --git a/src/Derivatives_2_1.F90 b/src/Derivatives_2_1.F90 index 9720dd6..ed97d9b 100644 --- a/src/Derivatives_2_1.F90 +++ b/src/Derivatives_2_1.F90 @@ -104,3 +104,186 @@ subroutine deriv_gf_2_1 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar ) end if end select direction end subroutine deriv_gf_2_1 + +subroutine up_deriv_gf_2_1 ( var, ni, nj, nk, dir, bb, gsize, delta, up, dvar ) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL, parameter :: zero = 0.0 + integer, parameter :: wp = kind(zero) + CCTK_INT, intent(IN) :: ni, nj, nk + CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var + CCTK_INT, intent(IN) :: dir + CCTK_INT, intent(IN) :: bb(2) + CCTK_INT, intent(IN) :: gsize + CCTK_REAL, intent(IN) :: delta + CCTK_INT, intent(IN) :: up + CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar + + CCTK_REAL, dimension(-1:1), save :: a1, a2 + CCTK_REAL, dimension(-1:1) :: a + CCTK_REAL, dimension(3,2), save :: q1, q2 + CCTK_REAL, dimension(3,2) :: q + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + a1(-1) = -1.000000000000000000000000000000000000000_wp + a1(0) = 1.000000000000000000000000000000000000000_wp + a1(1) = 0_wp + + q1(1,1) = 0_wp + q1(2,1) = 0_wp + q1(3,1) = 0_wp + q1(1,2) = -1.000000000000000000000000000000000000000_wp + q1(2,2) = 1.000000000000000000000000000000000000000_wp + q1(3,2) = 0_wp + + a2(-1) = 0_wp + a2(0) = -1.000000000000000000000000000000000000000_wp + a2(1) = 1.000000000000000000000000000000000000000_wp + + q2(1,1) = -2.000000000000000000000000000000000000000_wp + q2(2,1) = 2.000000000000000000000000000000000000000_wp + q2(3,1) = 0_wp + q2(1,2) = 0_wp + q2(2,2) = -1.000000000000000000000000000000000000000_wp + q2(3,2) = 1.000000000000000000000000000000000000000_wp + + first = .false. + end if + + idel = 1.0_wp / delta + + upwinding: select case (up) + case (-1) upwinding + a = a1 + case (1) upwinding + a = a2 + end select upwinding + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + upwinding_xl: select case (up) + case (-1) upwinding_xl + q = q1 + case (1) upwinding_xl + q = q2 + end select upwinding_xl + + dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) ) * idel + dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + & + q(3,2) * var(3,:,:) ) * idel + il = 3 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + upwinding_xr: select case (up) + case (-1) upwinding_xr + q = q2 + case (1) upwinding_xr + q = q1 + end select upwinding_xr + + dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + & + q(2,2) * var(ni-1,:,:) + & + q(3,2) * var(ni-2,:,:) ) * idel + dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + & + q(2,1) * var(ni-1,:,:) ) * idel + ir = ni - 2 + end if + dvar(il:ir,:,:) = ( a(-1) * var(il-1:ir-1,:,:) + & + a(0) * var(il:ir,:,:) + & + a(1) * var(il+1:ir+1,:,:) ) * idel + case (1) direction + if ( zero_derivs_y /= 0 ) then + dvar = zero + else + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + upwinding_yl: select case (up) + case (-1) upwinding_yl + q = q1 + case (1) upwinding_yl + q = q2 + end select upwinding_yl + + dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) ) * idel + dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + & + q(3,2) * var(:,3,:) ) * idel + jl = 3 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + upwinding_yr: select case (up) + case (-1) upwinding_yr + q = q2 + case (1) upwinding_yr + q = q1 + end select upwinding_yr + + dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + & + q(2,2) * var(:,nj-1,:) + & + q(3,2) * var(:,nj-2,:) ) * idel + dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + & + q(2,1) * var(:,nj-1,:) ) * idel + jr = nj - 2 + end if + dvar(:,jl:jr,:) = ( a(-1) * var(:,jl-1:jr-1,:) + & + a(0) * var(:,jl:jr,:) + & + a(1) * var(:,jl+1:jr+1,:) ) * idel + end if + case (2) direction + if ( zero_derivs_z /= 0 ) then + dvar = zero + else + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + upwinding_zl: select case (up) + case (-1) upwinding_zl + q = q1 + case (1) upwinding_zl + q = q2 + end select upwinding_zl + + dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) ) * idel + dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + & + q(3,2) * var(:,:,3) ) * idel + kl = 3 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + upwinding_zr: select case (up) + case (-1) upwinding_zr + q = q2 + case (1) upwinding_zr + q = q1 + end select upwinding_zr + + dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + & + q(2,2) * var(:,:,nk-1) + & + q(3,2) * var(:,:,nk-2) ) * idel + dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + & + q(2,1) * var(:,:,nk-1) ) * idel + kr = nk - 2 + end if + dvar(:,:,kl:kr) = ( a(-1) * var(:,:,kl-1:kr-1) + & + a(0) * var(:,:,kl:kr) + & + a(1) * var(:,:,kl+1:kr+1) ) * idel + end if + end select direction +end subroutine up_deriv_gf_2_1 diff --git a/src/Derivatives_4_2.F90 b/src/Derivatives_4_2.F90 index 7fc40b1..85e7ab3 100644 --- a/src/Derivatives_4_2.F90 +++ b/src/Derivatives_4_2.F90 @@ -166,3 +166,295 @@ subroutine deriv_gf_4_2 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar ) end if end select direction end subroutine deriv_gf_4_2 + +subroutine up_deriv_gf_4_2 ( var, ni, nj, nk, dir, bb, gsize, delta, up, dvar ) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL, parameter :: zero = 0.0 + integer, parameter :: wp = kind(zero) + CCTK_INT, intent(IN) :: ni, nj, nk + CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var + CCTK_INT, intent(IN) :: dir + CCTK_INT, intent(IN) :: bb(2) + CCTK_INT, intent(IN) :: gsize + CCTK_REAL, intent(IN) :: delta + CCTK_INT, intent(IN) :: up + CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar + + CCTK_REAL, dimension(-2:2), save :: a1, a2 + CCTK_REAL, dimension(-2:2) :: a + CCTK_REAL, dimension(6,4), save :: q1, q2 + CCTK_REAL, dimension(6,4) :: q + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + a1(-2) = 0.1666666666666666666666666666666666666667_wp + a1(-1) = -1.000000000000000000000000000000000000000_wp + a1(0) = 0.5000000000000000000000000000000000000000_wp + a1(1) = 0.3333333333333333333333333333333333333333_wp + a1(2) = 0_wp + + q1(1,1) = -1.176470588235294117647058823529411764706_wp + q1(2,1) = 1.264705882352941176470588235294117647059_wp + q1(3,1) = 0_wp + q1(4,1) = -0.08823529411764705882352941176470588235294_wp + q1(5,1) = 0_wp + q1(6,1) = 0_wp + q1(1,2) = -0.6355932203389830508474576271186440677966_wp + q1(2,2) = 0.3389830508474576271186440677966101694915_wp + q1(3,2) = 0.2288135593220338983050847457627118644068_wp + q1(4,2) = 0.06779661016949152542372881355932203389831_wp + q1(5,2) = 0_wp + q1(6,2) = 0_wp + q1(1,3) = 0.1860465116279069767441860465116279069767_wp + q1(2,3) = -1.058139534883720930232558139534883720930_wp + q1(3,3) = 0.5581395348837209302325581395348837209302_wp + q1(4,3) = 0.3139534883720930232558139534883720930233_wp + q1(5,3) = 0_wp + q1(6,3) = 0_wp + q1(1,4) = 0.03061224489795918367346938775510204081633_wp + q1(2,4) = 0.08163265306122448979591836734693877551020_wp + q1(3,4) = -0.9285714285714285714285714285714285714286_wp + q1(4,4) = 0.4897959183673469387755102040816326530612_wp + q1(5,4) = 0.3265306122448979591836734693877551020408_wp + q1(6,4) = 0_wp + + a2(-2) = 0_wp + a2(-1) = -0.3333333333333333333333333333333333333333_wp + a2(0) = -0.5000000000000000000000000000000000000000_wp + a2(1) = 1.000000000000000000000000000000000000000_wp + a2(2) = -0.1666666666666666666666666666666666666667_wp + + q2(1,1) = -1.647058823529411764705882352941176470588_wp + q2(2,1) = 2.205882352941176470588235294117647058824_wp + q2(3,1) = -0.4705882352941176470588235294117647058824_wp + q2(4,1) = -0.08823529411764705882352941176470588235294_wp + q2(5,1) = 0_wp + q2(6,1) = 0_wp + q2(1,2) = -0.3644067796610169491525423728813559322034_wp + q2(2,2) = -0.3389830508474576271186440677966101694915_wp + q2(3,2) = 0.7711864406779661016949152542372881355932_wp + q2(4,2) = -0.06779661016949152542372881355932203389831_wp + q2(5,2) = 0_wp + q2(6,2) = 0_wp + q2(1,3) = 0_wp + q2(2,3) = -0.3139534883720930232558139534883720930233_wp + q2(3,3) = -0.5581395348837209302325581395348837209302_wp + q2(4,3) = 1.058139534883720930232558139534883720930_wp + q2(5,3) = -0.1860465116279069767441860465116279069767_wp + q2(6,3) = 0_wp + q2(1,4) = 0.03061224489795918367346938775510204081633_wp + q2(2,4) = -0.08163265306122448979591836734693877551020_wp + q2(3,4) = -0.2755102040816326530612244897959183673469_wp + q2(4,4) = -0.4897959183673469387755102040816326530612_wp + q2(5,4) = 0.9795918367346938775510204081632653061224_wp + q2(6,4) = -0.1632653061224489795918367346938775510204_wp + + first = .false. + end if + + idel = 1.0_wp / delta + + upwinding: select case (up) + case (-1) upwinding + a = a1 + case (1) upwinding + a = a2 + end select upwinding + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + upwinding_xl: select case (up) + case (-1) upwinding_xl + q = q1 + case (1) upwinding_xl + q = q2 + end select upwinding_xl + + dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + & + q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) ) * idel + dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + & + q(3,2) * var(3,:,:) + q(4,2) * var(4,:,:) ) * idel + dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + & + q(3,3) * var(3,:,:) + q(4,3) * var(4,:,:) + & + q(5,3) * var(5,:,:) ) * idel + dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + & + q(3,4) * var(3,:,:) + q(4,4) * var(4,:,:) + & + q(5,4) * var(5,:,:) + q(6,4) * var(6,:,:) ) * idel + il = 5 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + upwinding_xr: select case (up) + case (-1) upwinding_xr + q = q2 + case (1) upwinding_xr + q = q1 + end select upwinding_xr + + dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + & + q(2,4) * var(ni-1,:,:) + & + q(3,4) * var(ni-2,:,:) + & + q(4,4) * var(ni-3,:,:) + & + q(5,4) * var(ni-4,:,:) + & + q(6,4) * var(ni-5,:,:) ) * idel + dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + & + q(2,3) * var(ni-1,:,:) + & + q(3,3) * var(ni-2,:,:) + & + q(4,3) * var(ni-3,:,:) + & + q(5,3) * var(ni-4,:,:) ) * idel + dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + & + q(2,2) * var(ni-1,:,:) + & + q(3,2) * var(ni-2,:,:) + & + q(4,2) * var(ni-3,:,:) ) * idel + dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + & + q(2,1) * var(ni-1,:,:) + & + q(3,1) * var(ni-2,:,:) + & + q(4,1) * var(ni-3,:,:) ) * idel + ir = ni - 4 + end if + dvar(il:ir,:,:) = ( a(-2) * var(il-2:ir-2,:,:) + & + a(-1) * var(il-1:ir-1,:,:) + & + a(0) * var(il:ir,:,:) + & + a(1) * var(il+1:ir+1,:,:) + & + a(2) * var(il+2:ir+2,:,:) ) * idel + case (1) direction + if ( zero_derivs_y /= 0 ) then + dvar = zero + else + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + upwinding_yl: select case (up) + case (-1) upwinding_yl + q = q1 + case (1) upwinding_yl + q = q2 + end select upwinding_yl + + dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + & + q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) ) * idel + dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + & + q(3,2) * var(:,3,:) + q(4,2) * var(:,4,:) ) * idel + dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + & + q(3,3) * var(:,3,:) + q(4,3) * var(:,4,:) + & + q(5,3) * var(:,5,:) ) * idel + dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + & + q(3,4) * var(:,3,:) + q(4,4) * var(:,4,:) + & + q(5,4) * var(:,5,:) + q(6,4) * var(:,6,:) ) * idel + jl = 5 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + upwinding_yr: select case (up) + case (-1) upwinding_yr + q = q2 + case (1) upwinding_yr + q = q1 + end select upwinding_yr + + dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + & + q(2,4) * var(:,nj-1,:) + & + q(3,4) * var(:,nj-2,:) + & + q(4,4) * var(:,nj-3,:) + & + q(5,4) * var(:,nj-4,:) + & + q(6,4) * var(:,nj-5,:) ) * idel + dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + & + q(2,3) * var(:,nj-1,:) + & + q(3,3) * var(:,nj-2,:) + & + q(4,3) * var(:,nj-3,:) + & + q(5,3) * var(:,nj-4,:) ) * idel + dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + & + q(2,2) * var(:,nj-1,:) + & + q(3,2) * var(:,nj-2,:) + & + q(4,2) * var(:,nj-3,:) ) * idel + dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + & + q(2,1) * var(:,nj-1,:) + & + q(3,1) * var(:,nj-2,:) + & + q(4,1) * var(:,nj-3,:) ) * idel + jr = nj - 4 + end if + dvar(:,jl:jr,:) = ( a(-2) * var(:,jl-2:jr-2,:) + & + a(-1) * var(:,jl-1:jr-1,:) + & + a(0) * var(:,jl:jr,:) + & + a(1) * var(:,jl+1:jr+1,:) + & + a(2) * var(:,jl+2:jr+2,:) ) * idel + end if + case (2) direction + if ( zero_derivs_z /= 0 ) then + dvar = zero + else + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + upwinding_zl: select case (up) + case (-1) upwinding_zl + q = q1 + case (1) upwinding_zl + q = q2 + end select upwinding_zl + + dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + & + q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) ) * idel + dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + & + q(3,2) * var(:,:,3) + q(4,2) * var(:,:,4) ) * idel + dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + & + q(3,3) * var(:,:,3) + q(4,3) * var(:,:,4) + & + q(5,3) * var(:,:,5) ) * idel + dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + & + q(3,4) * var(:,:,3) + q(4,4) * var(:,:,4) + & + q(5,4) * var(:,:,5) + q(6,4) * var(:,:,6) ) * idel + kl = 5 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + upwinding_zr: select case (up) + case (-1) upwinding_zr + q = q2 + case (1) upwinding_zr + q = q1 + end select upwinding_zr + + dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + & + q(2,4) * var(:,:,nk-1) + & + q(3,4) * var(:,:,nk-2) + & + q(4,4) * var(:,:,nk-3) + & + q(5,4) * var(:,:,nk-4) + & + q(6,4) * var(:,:,nk-5) ) * idel + dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + & + q(2,3) * var(:,:,nk-1) + & + q(3,3) * var(:,:,nk-2) + & + q(4,3) * var(:,:,nk-3) + & + q(5,3) * var(:,:,nk-4) ) * idel + dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + & + q(2,2) * var(:,:,nk-1) + & + q(3,2) * var(:,:,nk-2) + & + q(4,2) * var(:,:,nk-3) ) * idel + dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + & + q(2,1) * var(:,:,nk-1) + & + q(3,1) * var(:,:,nk-2) + & + q(4,1) * var(:,:,nk-3) ) * idel + kl = nk - 4 + end if + dvar(:,:,kl:kr) = ( a(-2) * var(:,:,kl-2:kr-2) + & + a(-1) * var(:,:,kl-1:kr-1) + & + a(0) * var(:,:,kl:kr) + & + a(1) * var(:,:,kl+1:kr+1) + & + a(2) * var(:,:,kl+2:kr+2) ) * idel + end if + end select direction +end subroutine up_deriv_gf_4_2 diff --git a/src/Derivatives_6_3_min_err_coeff.F90 b/src/Derivatives_6_3_min_err_coeff.F90 index 703eb4e..ae26440 100644 --- a/src/Derivatives_6_3_min_err_coeff.F90 +++ b/src/Derivatives_6_3_min_err_coeff.F90 @@ -287,3 +287,429 @@ subroutine deriv_gf_6_3_opt ( var, ni, nj, nk, dir, bb, gsize, delta, dvar ) end if end select direction end subroutine deriv_gf_6_3_opt + +subroutine up_deriv_gf_6_3_opt ( var, ni, nj, nk, dir, bb, gsize, delta, up, dvar ) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL, parameter :: zero = 0.0 + integer, parameter :: wp = kind(zero) + CCTK_INT, intent(IN) :: ni, nj, nk + CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var + CCTK_INT, intent(IN) :: dir + CCTK_INT, intent(IN) :: bb(2) + CCTK_INT, intent(IN) :: gsize + CCTK_REAL, intent(IN) :: delta + CCTK_INT, intent(IN) :: up + CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar + + CCTK_REAL, dimension(-3:3), save :: a1, a2 + CCTK_REAL, dimension(-3:3), save :: a + CCTK_REAL, dimension(9,6), save :: q1, q2 + CCTK_REAL, dimension(9,6), save :: q + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + a1(-3) = -0.03333333333333333333333333333333333333333_wp + a1(-2) = 0.2500000000000000000000000000000000000000_wp + a1(-1) = -1.000000000000000000000000000000000000000_wp + a1(0) = 0.3333333333333333333333333333333333333333_wp + a1(1) = 0.5000000000000000000000000000000000000000_wp + a1(2) = -0.05000000000000000000000000000000000000000_wp + a1(3) = 0_wp + + q1(1,1) = -1.529782401641145871492417026888416733827_wp + q1(2,1) = 1.875173434574214612019635448027523684006_wp + q1(3,1) = 0.01654806207923748080513710988516743214616_wp + q1(4,1) = -0.5028607772715413996090537489364770622056_wp + q1(5,1) = 0.1042956382142412661862395105494407610836_wp + q1(6,1) = 0.03662604404499391209045870736276191879693_wp + q1(7,1) = 0_wp + q1(8,1) = 0_wp + q1(9,1) = 0_wp + q1(1,2) = -0.4980311697078740570957463452947252270539_wp + q1(2,2) = 0.1198701406809289935902771997003246482977_wp + q1(3,2) = 0.1435237728855054645449658808725554143554_wp + q1(4,2) = 0.3305195343366502157809026121276381376866_wp + q1(5,2) = -0.08093510151416251407907993920061497636080_wp + q1(6,2) = -0.01494717668104810274131940820517799692506_wp + q1(7,2) = 0_wp + q1(8,2) = 0_wp + q1(9,2) = 0_wp + q1(1,3) = 0.1510193840162481602541749755506215314879_wp + q1(2,3) = -0.9553948881729209047544587102401343033294_wp + q1(3,3) = 0.5046108447067502766506824050165990409443_wp + q1(4,3) = 0.2083429547547929960329552908847933822643_wp + q1(5,3) = 0.1359641889806059957290021663880456685165_wp + q1(6,3) = -0.04454248428547652391235612759992531988360_wp + q1(7,3) = 0_wp + q1(8,3) = 0_wp + q1(9,3) = 0_wp + q1(1,4) = 0.1012044551031772450693035010120167087525_wp + q1(2,4) = -0.2092303756284921666519857323651163414844_wp + q1(3,4) = -0.5084563818511371174184254139930350549204_wp + q1(4,4) = 0.2687068482925919014741556260496361261429_wp + q1(5,4) = 0.3188546759663099905683904142802147346506_wp + q1(6,4) = 0.02892077811755014695856160501628382685880_wp + q1(7,4) = 0_wp + q1(8,4) = 0_wp + q1(9,4) = 0_wp + q1(1,5) = -0.03614399304268576976452921364705641609825_wp + q1(2,5) = 0.08686979490791345456798112347556020198327_wp + q1(3,5) = 0.1257842030401363832877174373294422223313_wp + q1(4,5) = -0.9822882337192980168734300444655759205262_wp + q1(5,5) = 0.3656214294782277516821124793703186492319_wp + q1(6,5) = 0.4950000137574403598524650898428590604628_wp + q1(7,5) = -0.05484321442173416275231687190554779738479_wp + q1(8,5) = 0_wp + q1(9,5) = 0_wp + q1(1,6) = -0.01141318406360863692889821914555232596651_wp + q1(2,6) = 0.02049729840293952857599941220163960606616_wp + q1(3,6) = -0.005307030684704263298464084855734448481896_wp + q1(4,6) = 0.1618715440442110368368458160128159110211_wp + q1(5,6) = -0.9382337284956231267046263227657132050941_wp + q1(6,6) = 0.3287596173603342389443163398095933882788_wp + q1(7,6) = 0.4931394260405013584164745097143900824182_wp + q1(8,6) = -0.04931394260405013584164745097143900824182_wp + q1(9,6) = 0_wp + + a2(-3) = 0_wp + a2(-2) = 0.05000000000000000000000000000000000000000_wp + a2(-1) = -0.5000000000000000000000000000000000000000_wp + a2(0) = -0.3333333333333333333333333333333333333333_wp + a2(1) = 1.000000000000000000000000000000000000000_wp + a2(2) = -0.2500000000000000000000000000000000000000_wp + a2(3) = 0.03333333333333333333333333333333333333333_wp + + q2(1,1) = -1.635284636237086966078100959777273060297_wp + q2(2,1) = 2.191680138362037895776687246694092663418_wp + q2(3,1) = -0.2999586417085858029519146887814015472662_wp + q2(4,1) = -0.3973585426756003050233698160476207357348_wp + q2(5,1) = 0.1042956382142412661862395105494407610836_wp + q2(6,1) = 0.03662604404499391209045870736276191879693_wp + q2(7,1) = 0_wp + q2(8,1) = 0_wp + q2(9,1) = 0_wp + q2(1,2) = -0.4261090852993166609415800254745304380753_wp + q2(2,2) = -0.1198701406809289935902771997003246482977_wp + q2(3,2) = 0.4312121105197350491616311601533345702699_wp + q2(4,2) = 0.1866753655195354234725699724872485597294_wp + q2(5,2) = -0.05696107337797671536102449926055004670126_wp + q2(6,2) = -0.01494717668104810274131940820517799692506_wp + q2(7,2) = 0_wp + q2(8,2) = 0_wp + q2(9,2) = 0_wp + q2(1,3) = -0.008331409049041400793408941823041323547141_wp + q2(2,3) = -0.3179917159117626605641230407454828831892_wp + q2(3,3) = -0.5046108447067502766506824050165990409443_wp + q2(4,3) = 1.005096920081240801270874877753107657440_wp + q2(5,3) = -0.1827373971499731263661656683592800415536_wp + q2(6,3) = 0.008574446736286663103505178191295631794744_wp + q2(7,3) = 0_wp + q2(8,3) = 0_wp + q2(9,3) = 0_wp + q2(1,4) = 0.1280751399324364352167190636169803213668_wp + q2(2,4) = -0.3704544846040473075364791079948980171701_wp + q2(3,4) = -0.1053961094122492652071919749185808657060_wp + q2(4,4) = -0.2687068482925919014741556260496361261429_wp + q2(5,4) = 0.7219149484051978427796238533546689238650_wp + q2(6,4) = -0.1323033308580049939259317706134978488270_wp + q2(7,4) = 0.02687068482925919014741556260496361261429_wp + q2(8,4) = 0_wp + q2(9,4) = 0_wp + q2(1,5) = -0.03614399304268576976452921364705641609825_wp + q2(2,5) = 0.1234319378557362297361923714125920669065_wp + q2(3,5) = -0.09358865464680026772155005029274896720786_wp + q2(4,5) = -0.4338560895019563893502613254100979466783_wp + q2(5,5) = -0.3656214294782277516821124793703186492319_wp + q2(6,5) = 1.043432157974781987375633808898337034311_wp + q2(7,5) = -0.2742160721086708137615843595277389869240_wp + q2(8,5) = 0.03656214294782277516821124793703186492319_wp + q2(9,5) = 0_wp + q2(1,6) = -0.01141318406360863692889821914555232596651_wp + q2(2,6) = 0.02049729840293952857599941220163960606616_wp + q2(3,6) = 0.02756893105132916059596754912522489034598_wp + q2(4,6) = -0.03538422637198950652974398787294012194614_wp + q2(5,6) = -0.4450943024551217682881518130513231226759_wp + q2(6,6) = -0.3287596173603342389443163398095933882788_wp + q2(7,6) = 0.9862788520810027168329490194287801648364_wp + q2(8,6) = -0.2465697130202506792082372548571950412091_wp + q2(9,6) = 0.03287596173603342389443163398095933882788_wp + + first = .false. + end if + + idel = 1.0_wp / delta + + upwinding: select case (up) + case (-1) upwinding + a = a1 + case (1) upwinding + a = a2 + end select upwinding + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + upwinding_xl: select case (up) + case (-1) upwinding_xl + q = q1 + case (1) upwinding_xl + q = q2 + end select upwinding_xl + + dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + & + q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) + & + q(5,1) * var(5,:,:) + q(6,1) * var(6,:,:) ) * idel + dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + & + q(3,2) * var(3,:,:) + q(4,2) * var(4,:,:) + & + q(5,2) * var(5,:,:) + q(6,2) * var(6,:,:) ) * idel + dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + & + q(3,3) * var(3,:,:) + q(4,3) * var(4,:,:) + & + q(5,3) * var(5,:,:) + q(6,3) * var(6,:,:) ) * idel + dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + & + q(3,4) * var(3,:,:) + q(4,4) * var(4,:,:) + & + q(5,4) * var(5,:,:) + q(6,4) * var(6,:,:) + & + q(7,4) * var(7,:,:) ) * idel + dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + & + q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + & + q(5,5) * var(5,:,:) + q(6,5) * var(6,:,:) + & + q(7,5) * var(7,:,:) + q(8,5) * var(8,:,:) ) * idel + dvar(6,:,:) = ( q(1,6) * var(1,:,:) + q(2,6) * var(2,:,:) + & + q(3,6) * var(3,:,:) + q(4,6) * var(4,:,:) + & + q(5,6) * var(5,:,:) + q(6,6) * var(6,:,:) + & + q(7,6) * var(7,:,:) + q(8,6) * var(8,:,:) + & + q(9,6) * var(9,:,:) ) * idel + + il = 7 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + upwinding_xr: select case (up) + case (-1) upwinding_xr + q = q2 + case (1) upwinding_xr + q = q1 + end select upwinding_xr + + dvar(ni-5,:,:) = - ( q(1,6) * var(ni,:,:) + q(2,6) * var(ni-1,:,:) + & + q(3,6) * var(ni-2,:,:) + q(4,6) * var(ni-3,:,:) + & + q(5,6) * var(ni-4,:,:) + q(6,6) * var(ni-5,:,:) + & + q(7,6) * var(ni-6,:,:) + q(8,6) * var(ni-7,:,:) + & + q(9,6) * var(ni-8,:,:) ) * idel + dvar(ni-4,:,:) = - ( q(1,5) * var(ni,:,:) + q(2,5) * var(ni-1,:,:) + & + q(3,5) * var(ni-2,:,:) + q(4,5) * var(ni-3,:,:) + & + q(5,5) * var(ni-4,:,:) + q(6,5) * var(ni-5,:,:) + & + q(7,5) * var(ni-6,:,:) + & + q(8,5) * var(ni-7,:,:) ) * idel + dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + & + q(3,4) * var(ni-2,:,:) + q(4,4) * var(ni-3,:,:) + & + q(5,4) * var(ni-4,:,:) + q(6,4) * var(ni-5,:,:) + & + q(7,4) * var(ni-6,:,:) ) * idel + dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + & + q(3,3) * var(ni-2,:,:) + q(4,3) * var(ni-3,:,:) + & + q(5,3) * var(ni-4,:,:) + & + q(6,3) * var(ni-5,:,:) ) * idel + dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(2,2) * var(ni-1,:,:) + & + q(3,2) * var(ni-2,:,:) + q(4,2) * var(ni-3,:,:) + & + q(5,2) * var(ni-4,:,:) + & + q(6,2) * var(ni-5,:,:) ) * idel + dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + q(2,1) * var(ni-1,:,:) + & + q(3,1) * var(ni-2,:,:) + q(4,1) * var(ni-3,:,:) + & + q(5,1) * var(ni-4,:,:) + & + q(6,1) * var(ni-5,:,:) ) * idel + + ir = ni - 6 + end if + dvar(il:ir,:,:) = ( a(-3) * var(il-3:ir-3,:,:) + & + a(-2) * var(il-2:ir-2,:,:) + & + a(-1) * var(il-1:ir-1,:,:) + & + a(0) * var(il:ir,:,:) + & + a(1) * var(il+1:ir+1,:,:) + & + a(2) * var(il+2:ir+2,:,:) + & + a(3) * var(il+3:ir+3,:,:) ) * idel + case (1) direction + if ( zero_derivs_y /= 0 ) then + dvar = zero + else + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + upwinding_yl: select case (up) + case (-1) upwinding_yl + q = q1 + case (1) upwinding_yl + q = q2 + end select upwinding_yl + + dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + & + q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) + & + q(5,1) * var(:,5,:) + q(6,1) * var(:,6,:) ) * idel + dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + & + q(3,2) * var(:,3,:) + q(4,2) * var(:,4,:) + & + q(5,2) * var(:,5,:) + q(6,2) * var(:,6,:) ) * idel + dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + & + q(3,3) * var(:,3,:) + q(4,3) * var(:,4,:) + & + q(5,3) * var(:,5,:) + q(6,3) * var(:,6,:) ) * idel + dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + & + q(3,4) * var(:,3,:) + q(4,4) * var(:,4,:) + & + q(5,4) * var(:,5,:) + q(6,4) * var(:,6,:) + & + q(7,4) * var(:,7,:) ) * idel + dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + & + q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + & + q(5,5) * var(:,5,:) + q(6,5) * var(:,6,:) + & + q(7,5) * var(:,7,:) + q(8,5) * var(:,8,:) ) * idel + dvar(:,6,:) = ( q(1,6) * var(:,1,:) + q(2,6) * var(:,2,:) + & + q(3,6) * var(:,3,:) + q(4,6) * var(:,4,:) + & + q(5,6) * var(:,5,:) + q(6,6) * var(:,6,:) + & + q(7,6) * var(:,7,:) + q(8,6) * var(:,8,:) + & + q(9,6) * var(:,9,:) ) * idel + + jl = 7 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + upwinding_yr: select case (up) + case (-1) upwinding_yr + q = q2 + case (1) upwinding_yr + q = q1 + end select upwinding_yr + + dvar(:,nj-5,:) = - ( q(1,6) * var(:,nj,:) + q(2,6) * var(:,nj-1,:) + & + q(3,6) * var(:,nj-2,:) + q(4,6) * var(:,nj-3,:) + & + q(5,6) * var(:,nj-4,:) + q(6,6) * var(:,nj-5,:) + & + q(7,6) * var(:,nj-6,:) + q(8,6) * var(:,nj-7,:) + & + q(9,6) * var(:,nj-8,:) ) * idel + dvar(:,nj-4,:) = - ( q(1,5) * var(:,nj,:) + q(2,5) * var(:,nj-1,:) + & + q(3,5) * var(:,nj-2,:) + q(4,5) * var(:,nj-3,:) + & + q(5,5) * var(:,nj-4,:) + q(6,5) * var(:,nj-5,:) + & + q(7,5) * var(:,nj-6,:) + & + q(8,5) * var(:,nj-7,:) ) * idel + dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + & + q(3,4) * var(:,nj-2,:) + q(4,4) * var(:,nj-3,:) + & + q(5,4) * var(:,nj-4,:) + q(6,4) * var(:,nj-5,:) + & + q(7,4) * var(:,nj-6,:) ) * idel + dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + & + q(3,3) * var(:,nj-2,:) + q(4,3) * var(:,nj-3,:) + & + q(5,3) * var(:,nj-4,:) + & + q(6,3) * var(:,nj-5,:) ) * idel + dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(2,2) * var(:,nj-1,:) + & + q(3,2) * var(:,nj-2,:) + q(4,2) * var(:,nj-3,:) + & + q(5,2) * var(:,nj-4,:) + & + q(6,2) * var(:,nj-5,:) ) * idel + dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + q(2,1) * var(:,nj-1,:) + & + q(3,1) * var(:,nj-2,:) + q(4,1) * var(:,nj-3,:) + & + q(5,1) * var(:,nj-4,:) + & + q(6,1) * var(:,nj-5,:) ) * idel + + jr = nj - 6 + end if + dvar(:,jl:jr,:) = ( a(-3) * var(:,jl-3:jr-3,:) + & + a(-2) * var(:,jl-2:jr-2,:) + & + a(-1) * var(:,jl-1:jr-1,:) + & + a(0) * var(:,jl:jr,:) + & + a(1) * var(:,jl+1:jr+1,:) + & + a(2) * var(:,jl+2:jr+2,:) + & + a(3) * var(:,jl+3:jr+3,:) ) * idel + end if + case (2) direction + if ( zero_derivs_z /= 0 ) then + dvar = zero + else + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + upwinding_zl: select case (up) + case (-1) upwinding_zl + q = q1 + case (1) upwinding_zl + q = q2 + end select upwinding_zl + + dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + & + q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) + & + q(5,1) * var(:,:,5) + q(6,1) * var(:,:,6) ) * idel + dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + & + q(3,2) * var(:,:,3) + q(4,2) * var(:,:,4) + & + q(5,2) * var(:,:,5) + q(6,2) * var(:,:,6) ) * idel + dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + & + q(3,3) * var(:,:,3) + q(4,3) * var(:,:,4) + & + q(5,3) * var(:,:,5) + q(6,3) * var(:,:,6) ) * idel + dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + & + q(3,4) * var(:,:,3) + q(4,4) * var(:,:,4) + & + q(5,4) * var(:,:,5) + q(6,4) * var(:,:,6) + & + q(7,4) * var(:,:,7) ) * idel + dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + & + q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + & + q(5,5) * var(:,:,5) + q(6,5) * var(:,:,6) + & + q(7,5) * var(:,:,7) + q(8,5) * var(:,:,8) ) * idel + dvar(:,:,6) = ( q(1,6) * var(:,:,1) + q(2,6) * var(:,:,2) + & + q(3,6) * var(:,:,3) + q(4,6) * var(:,:,4) + & + q(5,6) * var(:,:,5) + q(6,6) * var(:,:,6) + & + q(7,6) * var(:,:,7) + q(8,6) * var(:,:,8) + & + q(9,6) * var(:,:,9) ) * idel + + kl = 7 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + upwinding_zr: select case (up) + case (-1) upwinding_zr + q = q2 + case (1) upwinding_zr + q = q1 + end select upwinding_zr + + dvar(:,:,nk-5) = - ( q(1,6) * var(:,:,nk) + q(2,6) * var(:,:,nk-1) + & + q(3,6) * var(:,:,nk-2) + q(4,6) * var(:,:,nk-3) + & + q(5,6) * var(:,:,nk-4) + q(6,6) * var(:,:,nk-5) + & + q(7,6) * var(:,:,nk-6) + q(8,6) * var(:,:,nk-7) + & + q(9,6) * var(:,:,nk-8) ) * idel + dvar(:,:,nk-4) = - ( q(1,5) * var(:,:,nk) + q(2,5) * var(:,:,nk-1) + & + q(3,5) * var(:,:,nk-2) + q(4,5) * var(:,:,nk-3) + & + q(5,5) * var(:,:,nk-4) + q(6,5) * var(:,:,nk-5) + & + q(7,5) * var(:,:,nk-6) + & + q(8,5) * var(:,:,nk-7) ) * idel + dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + & + q(3,4) * var(:,:,nk-2) + q(4,4) * var(:,:,nk-3) + & + q(5,4) * var(:,:,nk-4) + q(6,4) * var(:,:,nk-5) + & + q(7,4) * var(:,:,nk-6) ) * idel + dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + & + q(3,3) * var(:,:,nk-2) + q(4,3) * var(:,:,nk-3) + & + q(5,3) * var(:,:,nk-4) + & + q(6,3) * var(:,:,nk-5) ) * idel + dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(2,2) * var(:,:,nk-1) + & + q(3,2) * var(:,:,nk-2) + q(4,2) * var(:,:,nk-3) + & + q(5,2) * var(:,:,nk-4) + & + q(6,2) * var(:,:,nk-5) ) * idel + dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + q(2,1) * var(:,:,nk-1) + & + q(3,1) * var(:,:,nk-2) + q(4,1) * var(:,:,nk-3) + & + q(5,1) * var(:,:,nk-4) + & + q(6,1) * var(:,:,nk-5) ) * idel + + kr = nk - 6 + end if + dvar(:,:,kl:kr) = ( a(-3) * var(:,:,kl-3:kr-3) + & + a(-2) * var(:,:,kl-2:kr-2) + & + a(-1) * var(:,:,kl-1:kr-1) + & + a(0) * var(:,:,kl:kr) + & + a(1) * var(:,:,kl+1:kr+1) + & + a(2) * var(:,:,kl+2:kr+2) + & + a(3) * var(:,:,kl+3:kr+3) ) * idel + + end if + end select direction +end subroutine up_deriv_gf_6_3_opt diff --git a/src/Derivatives_8_4_min_err_coeff.F90 b/src/Derivatives_8_4_min_err_coeff.F90 index 42eae82..63390ae 100644 --- a/src/Derivatives_8_4_min_err_coeff.F90 +++ b/src/Derivatives_8_4_min_err_coeff.F90 @@ -429,3 +429,622 @@ subroutine deriv_gf_8_4_opt ( var, ni, nj, nk, dir, bb, gsize, delta, dvar ) end if end select direction end subroutine deriv_gf_8_4_opt + +subroutine up_deriv_gf_8_4_opt ( var, ni, nj, nk, dir, bb, gsize, delta, up, dvar ) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL, parameter :: zero = 0.0 + integer, parameter :: wp = kind(zero) + CCTK_INT, intent(IN) :: ni, nj, nk + CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var + CCTK_INT, intent(IN) :: dir + CCTK_INT, intent(IN) :: bb(2) + CCTK_INT, intent(IN) :: gsize + CCTK_REAL, intent(IN) :: delta + CCTK_INT, intent(IN) :: up + CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar + + CCTK_REAL, dimension(-4:4), save :: a1, a2 + CCTK_REAL, dimension(-4:4), save :: a + CCTK_REAL, dimension(12,8), save :: q1, q2 + CCTK_REAL, dimension(12,8), save :: q + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr, j, k + + logical, save :: first = .true. + + if ( first ) then + a1(-4) = 0.007142857142857142857142857142857142857143_wp + a1(-3) = -0.06666666666666666666666666666666666666667_wp + a1(-2) = 0.3000000000000000000000000000000000000000_wp + a1(-1) = -1.000000000000000000000000000000000000000_wp + a1(0) = 0.2500000000000000000000000000000000000000_wp + a1(1) = 0.6000000000000000000000000000000000000000_wp + a1(2) = -0.1000000000000000000000000000000000000000_wp + a1(3) = 0.009523809523809523809523809523809523809524_wp + a1(4) = 0_wp + + q1(1,1) = -1.683432578685956376544499542432311020539_wp + q1(2,1) = 2.206194911574014373726323572790957192479_wp + q1(3,1) = 0.01048240012328238642042347285056590509052_wp + q1(4,1) = -0.7964177278055455600674641925035752355813_wp + q1(5,1) = 0.04375559918421074726330123748705042737925_wp + q1(6,1) = 0.3015760468426855559502835132688146245819_wp + q1(7,1) = -0.05430139927462465206917597295687512552949_wp + q1(8,1) = -0.02785725195806647467919208850462676788091_wp + q1(9,1) = 0_wp + q1(10,1) = 0_wp + q1(11,1) = 0_wp + q1(12,1) = 0_wp + q1(1,2) = -0.4451389901846271204230470438663315357999_wp + q1(2,2) = 0.03979384217746298000464189244717138065595_wp + q1(3,2) = 0.04117629055092890361160914968898310705568_wp + q1(4,2) = 0.5200531283081606472663885367960872642342_wp + q1(5,2) = -0.04132093916257534380117029102933732094903_wp + q1(6,2) = -0.1571076966553020594965523100929511511562_wp + q1(7,2) = 0.02929412303865500620649355002594213960776_wp + q1(8,2) = 0.01325024192729698663163651603043611635147_wp + q1(9,2) = 0_wp + q1(10,2) = 0_wp + q1(11,2) = 0_wp + q1(12,2) = 0_wp + q1(1,3) = 0.1544592597396873173299548977341719367844_wp + q1(2,3) = -1.020860666961436408391054234781341175626_wp + q1(3,3) = 0.7352247001589518809305910119936419247628_wp + q1(4,3) = -0.1443377614671085870438942853092495421459_wp + q1(5,3) = 0.4029099592852542668089969628980513126845_wp + q1(6,3) = -0.1256847607545176255327746274725689708913_wp + q1(7,3) = -0.01218362768224537680217842793138395523727_wp + q1(8,3) = 0.01047289768141453270035870286867846966855_wp + q1(9,3) = 0_wp + q1(10,3) = 0_wp + q1(11,3) = 0_wp + q1(12,3) = 0_wp + q1(1,4) = 0.1147228529148158318034738218073197976023_wp + q1(2,4) = -0.3538780254601914958115698874206531961307_wp + q1(3,4) = -0.1858994785495101740456777336107169678288_wp + q1(4,4) = 0.1370483808799853486477767415913728669524_wp + q1(5,4) = 0.03972858582710281564215102156521588367881_wp + q1(6,4) = 0.3165092134256211709442618752193670178495_wp + q1(7,4) = -0.05069213274750316093259254157360772183295_wp + q1(8,4) = -0.01753939629032033624782329757829768029059_wp + q1(9,4) = 0_wp + q1(10,4) = 0_wp + q1(11,4) = 0_wp + q1(12,4) = 0_wp + q1(1,5) = -0.01395722952056041916292538026746562805620_wp + q1(2,5) = 0.01429909551156929498323180030121632055909_wp + q1(3,5) = 0.2332633686362482363893930312775209255800_wp + q1(4,5) = -1.142300242635165935993053851544552454320_wp + q1(5,5) = 0.6057550738330768177827040634066377989296_wp + q1(6,5) = 0.2825179352075121763892365873707254541218_wp + q1(7,5) = 0.03248906901967492013975174489420235738435_wp + q1(8,5) = -0.01206707005235509052833799543828477419916_wp + q1(9,5) = 0_wp + q1(10,5) = 0_wp + q1(11,5) = 0_wp + q1(12,5) = 0_wp + q1(1,6) = -0.06956044896551484587496684695381321009857_wp + q1(2,6) = 0.1930764794969705526456194423505107491362_wp + q1(3,6) = -0.01938619631973248083178041998353003089641_wp + q1(4,6) = -0.2887164590533642022788579933885334942677_wp + q1(5,6) = -0.4040701147451219074780049288053335780831_wp + q1(6,6) = 0.1955440022503079624043376625896000853565_wp + q1(7,6) = 0.4088327186474708718613662440475307694516_wp + q1(8,6) = -0.02316927663483720615835459462174938908831_wp + q1(9,6) = 0.007449295323821255710641434765318098489772_wp + q1(10,6) = 0_wp + q1(11,6) = 0_wp + q1(12,6) = 0_wp + q1(1,7) = 0.01734328280101122760645584447264624051821_wp + q1(2,7) = -0.04840773491893482968357096782164797435370_wp + q1(3,7) = 0.01113356043999974094579726795398111697087_wp + q1(4,7) = 0.03683258325116868830333840353178374785556_wp + q1(5,7) = 0.2020929196538139300412937632456998385290_wp + q1(6,7) = -0.9993401510385675084229179036326626049945_wp + q1(7,7) = 0.2707691940348400441134114620253070105404_wp + q1(8,7) = 0.6075690064749917706800934212435751620821_wp + q1(9,7) = -0.1083076776139360176453645848101228042161_wp + q1(10,7) = 0.01031501691561295406146329379144026706820_wp + q1(11,7) = 0_wp + q1(12,7) = 0_wp + q1(1,8) = 0.008138876723484902470545040666116789759842_wp + q1(2,8) = -0.02002921745130282774392822052620050919550_wp + q1(3,8) = -0.002671343688148396036480262140129972226295_wp + q1(4,8) = 0.03832298452325717814060194729385919629895_wp + q1(5,8) = -0.05168030333817646678101782026786038350804_wp + q1(6,8) = 0.2274981280312773565363597229555297756596_wp + q1(7,8) = -0.9520792032950265164134838698189786312337_wp + q1(8,8) = 0.2476881136110866797237534905735256806167_wp + q1(9,8) = 0.5944514726666080313370083773764616334802_wp + q1(10,8) = -0.09907524544443467188950139622941027224670_wp + q1(11,8) = 0.009435737661374730656142990117086692594923_wp + q1(12,8) = 0_wp + + a2(-4) = 0_wp + a2(-3) = -0.009523809523809523809523809523809523809524_wp + a2(-2) = 0.1000000000000000000000000000000000000000_wp + a2(-1) = -0.6000000000000000000000000000000000000000_wp + a2(0) = -0.2500000000000000000000000000000000000000_wp + a2(1) = 1.000000000000000000000000000000000000000_wp + a2(2) = -0.3000000000000000000000000000000000000000_wp + a2(3) = 0.06666666666666666666666666666666666666667_wp + a2(4) = -0.007142857142857142857142857142857142857143_wp + + q2(1,1) = -1.707654630177840640955211766064430603569_wp + q2(2,1) = 2.303083117541551431369172467319435524597_wp + q2(3,1) = -0.1348499088280232000438498689421515930856_wp + q2(4,1) = -0.6995295218380085024246152979750969034639_wp + q2(5,1) = 0.01953354769232648285258901385493084434991_wp + q2(6,1) = 0.3015760468426855559502835132688146245819_wp + q2(7,1) = -0.05430139927462465206917597295687512552949_wp + q2(8,1) = -0.02785725195806647467919208850462676788091_wp + q2(9,1) = 0_wp + q2(10,1) = 0_wp + q2(11,1) = 0_wp + q2(12,1) = 0_wp + q2(1,2) = -0.4264124762187621886561567415382508860795_wp + q2(2,2) = -0.03979384217746298000464189244717138065595_wp + q2(3,2) = 0.1722618883119834259798412659855476550988_wp + q2(4,2) = 0.4170573014959035225484918739916436907718_wp + q2(5,2) = -0.003867911230845480267389686373176021508129_wp + q2(6,2) = -0.1617893251467682924382748856749713135863_wp + q2(7,2) = 0.02929412303865500620649355002594213960776_wp + q2(8,2) = 0.01325024192729698663163651603043611635147_wp + q2(9,2) = 0_wp + q2(10,2) = 0_wp + q2(11,2) = 0_wp + q2(12,2) = 0_wp + q2(1,3) = -0.01200671010762254250338646347193566882227_wp + q2(2,3) = -0.2440194743406570625021278824861723494614_wp + q2(3,3) = -0.7352247001589518809305910119936419247628_wp + q2(4,3) = 1.298367310542910198178397511810349706445_wp + q2(5,3) = -0.3739312333355250790799293893971175134799_wp + q2(6,3) = 0.09626986570856218757834718746890783658430_wp + q2(7,3) = -0.03992795599013035344106865479906855617172_wp + q2(8,3) = 0.01047289768141453270035870286867846966855_wp + q2(9,3) = 0_wp + q2(10,3) = 0_wp + q2(11,3) = 0_wp + q2(12,3) = 0_wp + q2(1,4) = 0.1306125202632199301974189512671891155098_wp + q2(2,4) = -0.4412711958764140369782680994499344446221_wp + q2(3,4) = 0.02066619697974310507560894936758416496903_wp + q2(4,4) = -0.1370483808799853486477767415913728669524_wp + q2(5,4) = 0.2621839287047601931573828340033863343841_wp + q2(6,4) = 0.2052815419867924821866459690002817924969_wp + q2(7,4) = -0.01891279805069496414470228265386908601791_wp + q2(8,4) = -0.02151181312742136084630957994326500976747_wp + q2(9,4) = 0_wp + q2(10,4) = 0_wp + q2(11,4) = 0_wp + q2(12,4) = 0_wp + q2(1,5) = -0.03126451734436261395671692493622670802562_wp + q2(2,5) = 0.1527573981019868533335641576513049603144_wp + q2(3,5) = -0.2513406904302132178367702194477893135637_wp + q2(4,5) = -0.1730921245022430275407273500939319760321_wp + q2(5,5) = -0.6057550738330768177827040634066377989296_wp + q2(6,5) = 1.251726053340435084841563088821345932409_wp + q2(7,5) = -0.4521149900467865340864115058311078817594_wp + q2(8,5) = 0.1263912325380624678219943619118038655562_wp + q2(9,5) = -0.01730728782380219479379154466876107996942_wp + q2(10,5) = 0_wp + q2(11,5) = 0_wp + q2(12,5) = 0_wp + q2(1,6) = -0.06956044896551484587496684695381321009857_wp + q2(2,6) = 0.1874895080041046108626383662765221752689_wp + q2(3,6) = 0.02530957562319505343206818860837856004222_wp + q2(4,6) = -0.4451516608536105722023281234602135625529_wp + q2(5,6) = -0.09119971114462916763106466866197344151272_wp + q2(6,6) = -0.1955440022503079624043376625896000853565_wp + q2(7,6) = 0.7217031222479636117083065041908909060220_wp + q2(8,6) = -0.1796044784350835760818247246934294573735_wp + q2(9,6) = 0.05214506726674878997449004335722668942840_wp + q2(10,6) = -0.005586971492865941782981076073988573867329_wp + q2(11,6) = 0_wp + q2(12,6) = 0_wp + q2(1,7) = 0.01734328280101122760645584447264624051821_wp + q2(2,7) = -0.04840773491893482968357096782164797435370_wp + q2(3,7) = 0.003397297753290025399699797610400916669720_wp + q2(4,7) = 0.09872268474484641267211816628042535026479_wp + q2(5,7) = -0.01452243557405810524943540637454576990329_wp + q2(6,7) = -0.5661094405828234378414595643921713881299_wp + q2(7,7) = -0.2707691940348400441134114620253070105404_wp + q2(8,7) = 1.040799716930735841261551760484066378947_wp + q2(9,7) = -0.3249230328418080529360937544303684126484_wp + q2(10,7) = 0.07220511840929067843024305654008186947743_wp + q2(11,7) = -0.007736262686709715546097470343580200301153_wp + q2(12,7) = 0_wp + q2(1,8) = 0.008138876723484902470545040666116789759842_wp + q2(2,8) = -0.02002921745130282774392822052620050919550_wp + q2(3,8) = -0.002671343688148396036480262140129972226295_wp + q2(4,8) = 0.03124618127722613014849470470604417685276_wp + q2(5,8) = 0.004934122630071917155840120434659772061498_wp + q2(6,8) = 0.02934763714240801275735693049670923116616_wp + q2(7,8) = -0.5557782215172878288554782849013375422469_wp + q2(8,8) = -0.2476881136110866797237534905735256806167_wp + q2(9,8) = 0.9907524544443467188950139622941027224670_wp + q2(10,8) = -0.2972257363333040156685041886882308167401_wp + q2(11,8) = 0.06605016362962311459300093081960684816446_wp + q2(12,8) = -0.007076803246031047992107242587815019446193_wp + + first = .false. + end if + + idel = 1.0_wp / delta + + upwinding: select case (up) + case (-1) upwinding + a = a1 + case (1) upwinding + a = a2 + end select upwinding + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + upwinding_xl: select case (up) + case (-1) upwinding_xl + q = q1 + case (1) upwinding_xl + q = q2 + end select upwinding_xl + + dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + & + q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) + & + q(5,1) * var(5,:,:) + q(6,1) * var(6,:,:) + & + q(7,1) * var(7,:,:) + q(8,1) * var(8,:,:) ) * idel + dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + & + q(3,2) * var(3,:,:) + q(4,2) * var(4,:,:) + & + q(5,2) * var(5,:,:) + q(6,2) * var(6,:,:) + & + q(7,2) * var(7,:,:) + q(8,2) * var(8,:,:) ) * idel + dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + & + q(3,3) * var(3,:,:) + q(4,3) * var(4,:,:) + & + q(5,3) * var(5,:,:) + q(6,3) * var(6,:,:) + & + q(7,3) * var(7,:,:) + q(8,3) * var(8,:,:) ) * idel + dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + & + q(3,4) * var(3,:,:) + q(4,4) * var(4,:,:) + & + q(5,4) * var(5,:,:) + q(6,4) * var(6,:,:) + & + q(7,4) * var(7,:,:) + q(8,4) * var(8,:,:) ) * idel + dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + & + q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + & + q(5,5) * var(5,:,:) + q(6,5) * var(6,:,:) + & + q(7,5) * var(7,:,:) + q(8,5) * var(8,:,:) + & + q(9,5) * var(9,:,:) ) * idel + dvar(6,:,:) = ( q(1,6) * var(1,:,:) + q(2,6) * var(2,:,:) + & + q(3,6) * var(3,:,:) + q(4,6) * var(4,:,:) + & + q(5,6) * var(5,:,:) + q(6,6) * var(6,:,:) + & + q(7,6) * var(7,:,:) + q(8,6) * var(8,:,:) + & + q(9,6) * var(9,:,:) + q(10,6) * var(10,:,:) ) * idel + dvar(7,:,:) = ( q(1,7) * var(1,:,:) + q(2,7) * var(2,:,:) + & + q(3,7) * var(3,:,:) + q(4,7) * var(4,:,:) + & + q(5,7) * var(5,:,:) + q(6,7) * var(6,:,:) + & + q(7,7) * var(7,:,:) + q(8,7) * var(8,:,:) + & + q(9,7) * var(9,:,:) + q(10,7) * var(10,:,:) + & + q(11,7) * var(11,:,:) ) * idel + dvar(8,:,:) = ( q(1,8) * var(1,:,:) + q(2,8) * var(2,:,:) + & + q(3,8) * var(3,:,:) + q(4,8) * var(4,:,:) + & + q(5,8) * var(5,:,:) + q(6,8) * var(6,:,:) + & + q(7,8) * var(7,:,:) + q(8,8) * var(8,:,:) + & + q(9,8) * var(9,:,:) + q(10,8) * var(10,:,:) + & + q(11,8) * var(11,:,:) + q(12,8) * var(12,:,:) ) * idel + il = 9 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + upwinding_xr: select case (up) + case (-1) upwinding_xr + q = q2 + case (1) upwinding_xr + q = q1 + end select upwinding_xr + + dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + q(2,1) * var(ni-1,:,:) + & + q(3,1) * var(ni-2,:,:) + q(4,1) * var(ni-3,:,:) + & + q(5,1) * var(ni-4,:,:) + q(6,1) * var(ni-5,:,:) + & + q(7,1) * var(ni-6,:,:) + & + q(8,1) * var(ni-7,:,:) ) * idel + dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(2,2) * var(ni-1,:,:) + & + q(3,2) * var(ni-2,:,:) + q(4,2) * var(ni-3,:,:) + & + q(5,2) * var(ni-4,:,:) + q(6,2) * var(ni-5,:,:) + & + q(7,2) * var(ni-6,:,:) + & + q(8,2) * var(ni-7,:,:) ) * idel + dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + & + q(3,3) * var(ni-2,:,:) + q(4,3) * var(ni-3,:,:) + & + q(5,3) * var(ni-4,:,:) + q(6,3) * var(ni-5,:,:) + & + q(7,3) * var(ni-6,:,:) + & + q(8,3) * var(ni-7,:,:) ) * idel + dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + & + q(3,4) * var(ni-2,:,:) + q(4,4) * var(ni-3,:,:) + & + q(5,4) * var(ni-4,:,:) + q(6,4) * var(ni-5,:,:) + & + q(7,4) * var(ni-6,:,:) + & + q(8,4) * var(ni-7,:,:) ) * idel + dvar(ni-4,:,:) = - ( q(1,5) * var(ni,:,:) + q(2,5) * var(ni-1,:,:) + & + q(3,5) * var(ni-2,:,:) + q(4,5) * var(ni-3,:,:) + & + q(5,5) * var(ni-4,:,:) + q(6,5) * var(ni-5,:,:) + & + q(7,5) * var(ni-6,:,:) + q(8,5) * var(ni-7,:,:) + & + q(9,5) * var(ni-8,:,:) ) * idel + dvar(ni-5,:,:) = - ( q(1,6) * var(ni,:,:) + q(2,6) * var(ni-1,:,:) + & + q(3,6) * var(ni-2,:,:) + q(4,6) * var(ni-3,:,:) + & + q(5,6) * var(ni-4,:,:) + q(6,6) * var(ni-5,:,:) + & + q(7,6) * var(ni-6,:,:) + q(8,6) * var(ni-7,:,:) + & + q(9,6) * var(ni-8,:,:) + & + q(10,6) * var(ni-9,:,:) ) * idel + dvar(ni-6,:,:) = - ( q(1,7) * var(ni,:,:) + q(2,7) * var(ni-1,:,:) + & + q(3,7) * var(ni-2,:,:) + q(4,7) * var(ni-3,:,:) + & + q(5,7) * var(ni-4,:,:) + q(6,7) * var(ni-5,:,:) + & + q(7,7) * var(ni-6,:,:) + q(8,7) * var(ni-7,:,:) + & + q(9,7) * var(ni-8,:,:) + q(10,7) * var(ni-9,:,:) + & + q(11,7) * var(ni-10,:,:) ) * idel + dvar(ni-7,:,:) = - ( q(1,8) * var(ni,:,:) + q(2,8) * var(ni-1,:,:) + & + q(3,8) * var(ni-2,:,:) + q(4,8) * var(ni-3,:,:) + & + q(5,8) * var(ni-4,:,:) + q(6,8) * var(ni-5,:,:) + & + q(7,8) * var(ni-6,:,:) + q(8,8) * var(ni-7,:,:) + & + q(9,8) * var(ni-8,:,:) + q(10,8) * var(ni-9,:,:) + & + q(11,8) * var(ni-10,:,:) + & + q(12,8) * var(ni-11,:,:) ) * idel + ir = ni - 8 + end if + dvar(il:ir,:,:) = ( a(-4) * var(il-4:ir-4,:,:) + & + a(-3) * var(il-3:ir-3,:,:) + & + a(-2) * var(il-2:ir-2,:,:) + & + a(-1) * var(il-1:ir-1,:,:) + & + a(0) * var(il:ir,:,:) + & + a(1) * var(il+1:ir+1,:,:) + & + a(2) * var(il+2:ir+2,:,:) + & + a(3) * var(il+3:ir+3,:,:) + & + a(4) * var(il+4:ir+4,:,:) ) * idel + case (1) direction + if ( zero_derivs_y /= 0 ) then + dvar = zero + else + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + upwinding_yl: select case (up) + case (-1) upwinding_yl + q = q1 + case (1) upwinding_yl + q = q2 + end select upwinding_yl + + dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + & + q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) + & + q(5,1) * var(:,5,:) + q(6,1) * var(:,6,:) + & + q(7,1) * var(:,7,:) + q(8,1) * var(:,8,:) ) * idel + dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + & + q(3,2) * var(:,3,:) + q(4,2) * var(:,4,:) + & + q(5,2) * var(:,5,:) + q(6,2) * var(:,6,:) + & + q(7,2) * var(:,7,:) + q(8,2) * var(:,8,:) ) * idel + dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + & + q(3,3) * var(:,3,:) + q(4,3) * var(:,4,:) + & + q(5,3) * var(:,5,:) + q(6,3) * var(:,6,:) + & + q(7,3) * var(:,7,:) + q(8,3) * var(:,8,:) ) * idel + dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + & + q(3,4) * var(:,3,:) + q(4,4) * var(:,4,:) + & + q(5,4) * var(:,5,:) + q(6,4) * var(:,6,:) + & + q(7,4) * var(:,7,:) + q(8,4) * var(:,8,:) ) * idel + dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + & + q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + & + q(5,5) * var(:,5,:) + q(6,5) * var(:,6,:) + & + q(7,5) * var(:,7,:) + q(8,5) * var(:,8,:) + & + q(9,5) * var(:,9,:) ) * idel + dvar(:,6,:) = ( q(1,6) * var(:,1,:) + q(2,6) * var(:,2,:) + & + q(3,6) * var(:,3,:) + q(4,6) * var(:,4,:) + & + q(5,6) * var(:,5,:) + q(6,6) * var(:,6,:) + & + q(7,6) * var(:,7,:) + q(8,6) * var(:,8,:) + & + q(9,6) * var(:,9,:) + q(10,6) * var(:,10,:) ) * idel + dvar(:,7,:) = ( q(1,7) * var(:,1,:) + q(2,7) * var(:,2,:) + & + q(3,7) * var(:,3,:) + q(4,7) * var(:,4,:) + & + q(5,7) * var(:,5,:) + q(6,7) * var(:,6,:) + & + q(7,7) * var(:,7,:) + q(8,7) * var(:,8,:) + & + q(9,7) * var(:,9,:) + q(10,7) * var(:,10,:) + & + q(11,7) * var(:,11,:) ) * idel + dvar(:,8,:) = ( q(1,8) * var(:,1,:) + q(2,8) * var(:,2,:) + & + q(3,8) * var(:,3,:) + q(4,8) * var(:,4,:) + & + q(5,8) * var(:,5,:) + q(6,8) * var(:,6,:) + & + q(7,8) * var(:,7,:) + q(8,8) * var(:,8,:) + & + q(9,8) * var(:,9,:) + q(10,8) * var(:,10,:) + & + q(11,8) * var(:,11,:) + q(12,8) * var(:,12,:) ) * idel + jl = 9 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + upwinding_yr: select case (up) + case (-1) upwinding_yr + q = q2 + case (1) upwinding_yr + q = q1 + end select upwinding_yr + + dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + q(2,1) * var(:,nj-1,:) + & + q(3,1) * var(:,nj-2,:) + q(4,1) * var(:,nj-3,:) + & + q(5,1) * var(:,nj-4,:) + q(6,1) * var(:,nj-5,:) + & + q(7,1) * var(:,nj-6,:) + & + q(8,1) * var(:,nj-7,:) ) * idel + dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(2,2) * var(:,nj-1,:) + & + q(3,2) * var(:,nj-2,:) + q(4,2) * var(:,nj-3,:) + & + q(5,2) * var(:,nj-4,:) + q(6,2) * var(:,nj-5,:) + & + q(7,2) * var(:,nj-6,:) + & + q(8,2) * var(:,nj-7,:) ) * idel + dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + & + q(3,3) * var(:,nj-2,:) + q(4,3) * var(:,nj-3,:) + & + q(5,3) * var(:,nj-4,:) + q(6,3) * var(:,nj-5,:) + & + q(7,3) * var(:,nj-6,:) + & + q(8,3) * var(:,nj-7,:) ) * idel + dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + & + q(3,4) * var(:,nj-2,:) + q(4,4) * var(:,nj-3,:) + & + q(5,4) * var(:,nj-4,:) + q(6,4) * var(:,nj-5,:) + & + q(7,4) * var(:,nj-6,:) + & + q(8,4) * var(:,nj-7,:) ) * idel + dvar(:,nj-4,:) = - ( q(1,5) * var(:,nj,:) + q(2,5) * var(:,nj-1,:) + & + q(3,5) * var(:,nj-2,:) + q(4,5) * var(:,nj-3,:) + & + q(5,5) * var(:,nj-4,:) + q(6,5) * var(:,nj-5,:) + & + q(7,5) * var(:,nj-6,:) + q(8,5) * var(:,nj-7,:) + & + q(9,5) * var(:,nj-8,:) ) * idel + dvar(:,nj-5,:) = - ( q(1,6) * var(:,nj,:) + q(2,6) * var(:,nj-1,:) + & + q(3,6) * var(:,nj-2,:) + q(4,6) * var(:,nj-3,:) + & + q(5,6) * var(:,nj-4,:) + q(6,6) * var(:,nj-5,:) + & + q(7,6) * var(:,nj-6,:) + q(8,6) * var(:,nj-7,:) + & + q(9,6) * var(:,nj-8,:) + & + q(10,6) * var(:,nj-9,:) ) * idel + dvar(:,nj-6,:) = - ( q(1,7) * var(:,nj,:) + q(2,7) * var(:,nj-1,:) + & + q(3,7) * var(:,nj-2,:) + q(4,7) * var(:,nj-3,:) + & + q(5,7) * var(:,nj-4,:) + q(6,7) * var(:,nj-5,:) + & + q(7,7) * var(:,nj-6,:) + q(8,7) * var(:,nj-7,:) + & + q(9,7) * var(:,nj-8,:) + & + q(10,7) * var(:,nj-9,:) + & + q(11,7) * var(:,nj-10,:) ) * idel + dvar(:,nj-7,:) = - ( q(1,8) * var(:,nj,:) + q(2,8) * var(:,nj-1,:) + & + q(3,8) * var(:,nj-2,:) + q(4,8) * var(:,nj-3,:) + & + q(5,8) * var(:,nj-4,:) + q(6,8) * var(:,nj-5,:) + & + q(7,8) * var(:,nj-6,:) + q(8,8) * var(:,nj-7,:) + & + q(9,8) * var(:,nj-8,:) + & + q(10,8) * var(:,nj-9,:) + & + q(11,8) * var(:,nj-10,:) + & + q(12,8) * var(:,nj-11,:) ) * idel + jr = nj - 8 + end if + dvar(:,jl:jr,:) = ( a(-4) * var(:,jl-4:jr-4,:) + & + a(-3) * var(:,jl-3:jr-3,:) + & + a(-2) * var(:,jl-2:jr-2,:) + & + a(-1) * var(:,jl-1:jr-1,:) + & + a(0) * var(:,jl:jr,:) + & + a(1) * var(:,jl+1:jr+1,:) + & + a(2) * var(:,jl+2:jr+2,:) + & + a(3) * var(:,jl+3:jr+3,:) + & + a(4) * var(:,jl+4:jr+4,:) ) * idel + end if + case (2) direction + if ( zero_derivs_z /= 0 ) then + dvar = zero + else + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + upwinding_zl: select case (up) + case (-1) upwinding_zl + q = q1 + case (1) upwinding_zl + q = q2 + end select upwinding_zl + + dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + & + q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) + & + q(5,1) * var(:,:,5) + q(6,1) * var(:,:,6) + & + q(7,1) * var(:,:,7) + q(8,1) * var(:,:,8) ) * idel + dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + & + q(3,2) * var(:,:,3) + q(4,2) * var(:,:,4) + & + q(5,2) * var(:,:,5) + q(6,2) * var(:,:,6) + & + q(7,2) * var(:,:,7) + q(8,2) * var(:,:,8) ) * idel + dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + & + q(3,3) * var(:,:,3) + q(4,3) * var(:,:,4) + & + q(5,3) * var(:,:,5) + q(6,3) * var(:,:,6) + & + q(7,3) * var(:,:,7) + q(8,3) * var(:,:,8) ) * idel + dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + & + q(3,4) * var(:,:,3) + q(4,4) * var(:,:,4) + & + q(5,4) * var(:,:,5) + q(6,4) * var(:,:,6) + & + q(7,4) * var(:,:,7) + q(8,4) * var(:,:,8) ) * idel + dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + & + q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + & + q(5,5) * var(:,:,5) + q(6,5) * var(:,:,6) + & + q(7,5) * var(:,:,7) + q(8,5) * var(:,:,8) + & + q(9,5) * var(:,:,9) ) * idel + dvar(:,:,6) = ( q(1,6) * var(:,:,1) + q(2,6) * var(:,:,2) + & + q(3,6) * var(:,:,3) + q(4,6) * var(:,:,4) + & + q(5,6) * var(:,:,5) + q(6,6) * var(:,:,6) + & + q(7,6) * var(:,:,7) + q(8,6) * var(:,:,8) + & + q(9,6) * var(:,:,9) + q(10,6) * var(:,:,10) ) * idel + dvar(:,:,7) = ( q(1,7) * var(:,:,1) + q(2,7) * var(:,:,2) + & + q(3,7) * var(:,:,3) + q(4,7) * var(:,:,4) + & + q(5,7) * var(:,:,5) + q(6,7) * var(:,:,6) + & + q(7,7) * var(:,:,7) + q(8,7) * var(:,:,8) + & + q(9,7) * var(:,:,9) + q(10,7) * var(:,:,10) + & + q(11,7) * var(:,:,11) ) * idel + dvar(:,:,8) = ( q(1,8) * var(:,:,1) + q(2,8) * var(:,:,2) + & + q(3,8) * var(:,:,3) + q(4,8) * var(:,:,4) + & + q(5,8) * var(:,:,5) + q(6,8) * var(:,:,6) + & + q(7,8) * var(:,:,7) + q(8,8) * var(:,:,8) + & + q(9,8) * var(:,:,9) + q(10,8) * var(:,:,10) + & + q(11,8) * var(:,:,11) + q(12,8) * var(:,:,12) ) * idel + kl = 9 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + upwinding_zr: select case (up) + case (-1) upwinding_zr + q = q2 + case (1) upwinding_zr + q = q1 + end select upwinding_zr + + dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + q(2,1) * var(:,:,nk-1) + & + q(3,1) * var(:,:,nk-2) + q(4,1) * var(:,:,nk-3) + & + q(5,1) * var(:,:,nk-4) + q(6,1) * var(:,:,nk-5) + & + q(7,1) * var(:,:,nk-6) + & + q(8,1) * var(:,:,nk-7) ) * idel + dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(2,2) * var(:,:,nk-1) + & + q(3,2) * var(:,:,nk-2) + q(4,2) * var(:,:,nk-3) + & + q(5,2) * var(:,:,nk-4) + q(6,2) * var(:,:,nk-5) + & + q(7,2) * var(:,:,nk-6) + & + q(8,2) * var(:,:,nk-7) ) * idel + dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + & + q(3,3) * var(:,:,nk-2) + q(4,3) * var(:,:,nk-3) + & + q(5,3) * var(:,:,nk-4) + q(6,3) * var(:,:,nk-5) + & + q(7,3) * var(:,:,nk-6) + & + q(8,3) * var(:,:,nk-7) ) * idel + dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + & + q(3,4) * var(:,:,nk-2) + q(4,4) * var(:,:,nk-3) + & + q(5,4) * var(:,:,nk-4) + q(6,4) * var(:,:,nk-5) + & + q(7,4) * var(:,:,nk-6) + & + q(8,4) * var(:,:,nk-7) ) * idel + dvar(:,:,nk-4) = - ( q(1,5) * var(:,:,nk) + q(2,5) * var(:,:,nk-1) + & + q(3,5) * var(:,:,nk-2) + q(4,5) * var(:,:,nk-3) + & + q(5,5) * var(:,:,nk-4) + q(6,5) * var(:,:,nk-5) + & + q(7,5) * var(:,:,nk-6) + q(8,5) * var(:,:,nk-7) + & + q(9,5) * var(:,:,nk-8) ) * idel + dvar(:,:,nk-5) = - ( q(1,6) * var(:,:,nk) + q(2,6) * var(:,:,nk-1) + & + q(3,6) * var(:,:,nk-2) + q(4,6) * var(:,:,nk-3) + & + q(5,6) * var(:,:,nk-4) + q(6,6) * var(:,:,nk-5) + & + q(7,6) * var(:,:,nk-6) + q(8,6) * var(:,:,nk-7) + & + q(9,6) * var(:,:,nk-8) + & + q(10,6) * var(:,:,nk-9) ) * idel + dvar(:,:,nk-6) = - ( q(1,7) * var(:,:,nk) + q(2,7) * var(:,:,nk-1) + & + q(3,7) * var(:,:,nk-2) + q(4,7) * var(:,:,nk-3) + & + q(5,7) * var(:,:,nk-4) + q(6,7) * var(:,:,nk-5) + & + q(7,7) * var(:,:,nk-6) + q(8,7) * var(:,:,nk-7) + & + q(9,7) * var(:,:,nk-8) + & + q(10,7) * var(:,:,nk-9) + & + q(11,7) * var(:,:,nk-10) ) * idel + dvar(:,:,nk-7) = - ( q(1,8) * var(:,:,nk) + q(2,8) * var(:,:,nk-1) + & + q(3,8) * var(:,:,nk-2) + q(4,8) * var(:,:,nk-3) + & + q(5,8) * var(:,:,nk-4) + q(6,8) * var(:,:,nk-5) + & + q(7,8) * var(:,:,nk-6) + q(8,8) * var(:,:,nk-7) + & + q(9,8) * var(:,:,nk-8) + & + q(10,8) * var(:,:,nk-9) + & + q(11,8) * var(:,:,nk-10) + & + q(12,8) * var(:,:,nk-11) ) * idel + kr = nk - 8 + end if + dvar(:,:,kl:kr) = ( a(-4) * var(:,:,kl-4:kr-4) + & + a(-3) * var(:,:,kl-3:kr-3) + & + a(-2) * var(:,:,kl-2:kr-2) + & + a(-1) * var(:,:,kl-1:kr-1) + & + a(0) * var(:,:,kl:kr) + & + a(1) * var(:,:,kl+1:kr+1) + & + a(2) * var(:,:,kl+2:kr+2) + & + a(3) * var(:,:,kl+3:kr+3) + & + a(4) * var(:,:,kl+4:kr+4) ) * idel + end if + end select direction +end subroutine up_deriv_gf_8_4_opt diff --git a/src/call_up_derivs.c b/src/call_up_derivs.c new file mode 100644 index 0000000..2189e3f --- /dev/null +++ b/src/call_up_derivs.c @@ -0,0 +1,206 @@ +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +void DiffUpGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, + const CCTK_INT up, const CCTK_REAL *var, CCTK_REAL *dvar, + const CCTK_INT table_handle ) +{ + cGH const * restrict const cctkGH = cctkGH_; + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + + CCTK_INT ni, nj, nk, gsize, ic, loc_order; + CCTK_REAL delta; + CCTK_INT ierr; + CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3]; + CCTK_INT symtable, pen_sym_handle; + CCTK_INT symbnd[6]; + + int nelements; + void CCTK_FCALL CCTK_FNAME(up_deriv_gf_2_1)(const CCTK_REAL *var, + const CCTK_INT *ni, + const CCTK_INT *nj, + const CCTK_INT *nk, + const CCTK_INT *dir, + const CCTK_INT *bb, + const CCTK_INT *gsize, + const CCTK_REAL *delta, + const CCTK_INT *up, + CCTK_REAL *dvar); + void CCTK_FCALL CCTK_FNAME(up_deriv_gf_4_2)(const CCTK_REAL *var, + const CCTK_INT *ni, + const CCTK_INT *nj, + const CCTK_INT *nk, + const CCTK_INT *dir, + const CCTK_INT *bb, + const CCTK_INT *gsize, + const CCTK_REAL *delta, + const CCTK_INT *up, + CCTK_REAL *dvar); + void CCTK_FCALL CCTK_FNAME(up_deriv_gf_6_3_opt)(const CCTK_REAL *var, + const CCTK_INT *ni, + const CCTK_INT *nj, + const CCTK_INT *nk, + const CCTK_INT *dir, + const CCTK_INT *bb, + const CCTK_INT *gsize, + const CCTK_REAL *delta, + const CCTK_INT *up, + CCTK_REAL *dvar); + void CCTK_FCALL CCTK_FNAME(up_deriv_gf_8_4_opt)(const CCTK_REAL *var, + const CCTK_INT *ni, + const CCTK_INT *nj, + const CCTK_INT *nk, + const CCTK_INT *dir, + const CCTK_INT *bb, + const CCTK_INT *gsize, + const CCTK_REAL *delta, + const CCTK_INT *up, + CCTK_REAL *dvar); + + + ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2]; + + if ( table_handle >=0 ) { + nelements = Util_TableGetInt ( table_handle, &loc_order, "order" ); + if ( nelements == UTIL_ERROR_TABLE_NO_SUCH_KEY ) { + loc_order = order; + } else if ( nelements != 1) { + CCTK_WARN (0, "The options table has an entry \"order\", but it does not have the right properties"); + } + } else { + loc_order = order; + } + + symtable = SymmetryTableHandleForGrid (cctkGH); + if (symtable<0) { + CCTK_WARN(0,"symtable is out of bounds"); + } + + ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle"); + if (ierr!=6) { + CCTK_WARN(0,"Util_TableGetIntArray returned error"); + } + + pen_sym_handle = SymmetryHandleOfName ( "multipatch" ); + + switch(dir) { + case 0: { + delta = CCTK_DELTA_SPACE(0); + if ( symbnd[0] == pen_sym_handle && cctk_bbox[0] == 1 ) { + bb[0] = 1; + } else { + bb[0] = 0; + } + if ( symbnd[1] == pen_sym_handle && cctk_bbox[1] == 1 ) { + bb[1] = 1; + } else { + bb[1] = 0; + } + gsize = cctk_nghostzones[0]; + break; + } + case 1: { + delta = CCTK_DELTA_SPACE(1); + if ( symbnd[2] == pen_sym_handle && cctk_bbox[2] == 1) { + bb[0] = 1; + } else { + bb[0] = 0; + } + if ( symbnd[3] == pen_sym_handle && cctk_bbox[3] == 1) { + bb[1] = 1; + } else { + bb[1] = 0; + } + gsize = cctk_nghostzones[1]; + break; + } + case 2: { + delta = CCTK_DELTA_SPACE(2); + if ( symbnd[4] == pen_sym_handle && cctk_bbox[4] == 1) { + bb[0] = 1; + } else { + bb[0] = 0; + } + if ( symbnd[4] == pen_sym_handle && cctk_bbox[5] == 1) { + bb[1] = 1; + } else { + bb[1] = 0; + } + gsize = cctk_nghostzones[2]; + break; + } + default: + CCTK_WARN (0, "Wrong direction specified"); + } + + if ( CCTK_Equals(norm_type,"Diagonal") ) { + if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) { + switch(loc_order) { + case 2: { + CCTK_FNAME(up_deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,&up,dvar); + break; + } + case 4: { + CCTK_FNAME(up_deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,&up,dvar); + break; + } + default: + CCTK_WARN (0, "Unknown stencil specified"); + } + } else { + switch(loc_order) { + case 2: { + CCTK_FNAME(up_deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,&up,dvar); + break; + } + case 4: { + CCTK_FNAME(up_deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,&up,dvar); + break; + } + case 6: { + CCTK_FNAME(up_deriv_gf_6_3_opt)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,&up,dvar); + break; + } + case 8: { + CCTK_FNAME(up_deriv_gf_8_4_opt)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,&up,dvar); + break; + } + default: + CCTK_WARN (0, "Unknown stencil specified"); + } + } + } else { + if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) { + switch(loc_order) { +/* case 4: { + CCTK_FNAME(deriv_gf_4_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar); + break; + } + case 6: { + CCTK_FNAME(deriv_gf_6_5)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar); + break; + } */ + default: + CCTK_WARN (0, "Unknown stencil specified"); + } + } else { + switch(loc_order) { +/* case 4: { + CCTK_FNAME(deriv_gf_4_3_opt)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar); + break; + } + case 6: { + CCTK_FNAME(deriv_gf_6_5_opt)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar); + break; + } */ + default: + CCTK_WARN (0, "Unknown stencil specified"); + } + } + } +} diff --git a/src/make.code.defn b/src/make.code.defn index 0371d01..dd2255b 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -3,6 +3,7 @@ # Source files in this directory SRCS = call_derivs.c \ + call_up_derivs.c \ call_derivs_name.c \ Derivatives_2_1.F90 \ Derivatives_4_2.F90 \ -- cgit v1.2.3