From b613d73f118c88bd116f0d8fe1c32cd3a7f855db Mon Sep 17 00:00:00 2001 From: diener Date: Thu, 8 May 2008 20:28:05 +0000 Subject: Add non-sbp upwinding operators. These have the same order as the centered operators but does not satisfy any SBP property. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@114 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- param.ccl | 4 + src/All_Coeffs_mod.F90 | 273 +++++++++++++++++++++++++++++++++ src/Coefficients_2_1.F90 | 82 ++++++++++ src/Coefficients_4_2.F90 | 94 ++++++++++++ src/Coefficients_6_3_min_err_coeff.F90 | 106 +++++++++++++ src/Coefficients_8_4_min_err_coeff.F90 | 118 ++++++++++++++ src/get_up_coeffs.c | 90 +++++++++-- 7 files changed, 750 insertions(+), 17 deletions(-) diff --git a/param.ccl b/param.ccl index 2ffd8d2..eba0e40 100644 --- a/param.ccl +++ b/param.ccl @@ -17,6 +17,10 @@ BOOLEAN sbp_2nd_deriv "Should the 2nd derivative operator be SBP" STEERABLE=alwa { } "yes" +BOOLEAN sbp_upwind_deriv "Should the upwind derivative operator be SBP" STEERABLE=always +{ +} "yes" + INT order "Order of accuracy" STEERABLE=always { 2:8:2 :: "" diff --git a/src/All_Coeffs_mod.F90 b/src/All_Coeffs_mod.F90 index 539e46e..1886497 100644 --- a/src/All_Coeffs_mod.F90 +++ b/src/All_Coeffs_mod.F90 @@ -59,6 +59,41 @@ module All_Coeffs_mod end subroutine coeffs_up_2_1 + subroutine coeffs_up_2 ( a1, q1, a2, q2 ) + + CCTK_REAL, dimension(-2:2), intent(OUT) :: a1, a2 + CCTK_REAL, dimension(3,2), intent(OUT) :: q1, q2 + + a1(-2) = 0.5_wp + a1(-1) = -2.0_wp + a1(0) = 1.5_wp + a1(1) = zero + a1(2) = zero + + q1(1,1) = zero + q1(2,1) = zero + q1(3,1) = zero + + q1(1,2) = -1.0_wp + q1(2,2) = 1.0_wp + q1(3,2) = zero + + a2(-2) = zero + a2(-1) = zero + a2(0) = -1.5_wp + a2(1) = 2.0_wp + a2(2) = -0.5_wp + + q2(1,1) = zero + q2(2,1) = zero + q2(3,1) = zero + + q2(1,2) = zero + q2(2,2) = -1.0_wp + q2(3,2) = 1.0_wp + + end subroutine coeffs_up_2 + subroutine coeffs_2_2_1 ( a, q ) CCTK_REAL, dimension(2), intent(OUT) :: a @@ -201,6 +236,59 @@ module All_Coeffs_mod end subroutine coeffs_up_4_2 + subroutine coeffs_up_4 ( a1, q1, a2, q2 ) + + CCTK_REAL, dimension(-3:3), intent(OUT) :: a1, a2 + CCTK_REAL, dimension(4,3), intent(OUT) :: q1, q2 + + a1(-3) = -0.08333333333333333333333333333333333333333_wp + a1(-2) = 0.5_wp + a1(-1) = -1.5_wp + a1(0) = 0.8333333333333333333333333333333333333333_wp + a1(1) = 0.25_wp + a1(2) = zero + a1(3) = zero + + q1(1,1) = zero + q1(2,1) = zero + q1(3,1) = zero + q1(4,1) = zero + + q1(1,2) = -1.0_wp + q1(2,2) = 1.0_wp + q1(3,2) = zero + q1(4,2) = zero + + q1(1,3) = 0.5_wp + q1(2,3) = -2.0_wp + q1(3,3) = 1.5_wp + q1(4,3) = zero + + a2(-3) = zero + a2(-2) = zero + a2(-1) = -0.25_wp + a2(0) = -0.8333333333333333333333333333333333333333_wp + a2(1) = 1.5_wp + a2(2) = -0.5_wp + a2(3) = 0.08333333333333333333333333333333333333333_wp + + q2(1,1) = zero + q2(2,1) = zero + q2(3,1) = zero + q2(4,1) = zero + + q2(1,2) = zero + q2(2,2) = zero + q2(3,2) = -1.0_wp + q2(4,2) = 1.0_wp + + q2(1,3) = zero + q2(2,3) = -1.5_wp + q2(3,3) = 2.0_wp + q2(4,3) = -0.5_wp + + end subroutine coeffs_up_4 + subroutine coeffs_2_4_2 ( a, q ) CCTK_REAL, dimension(3), intent(OUT) :: a @@ -581,6 +669,81 @@ module All_Coeffs_mod end subroutine coeffs_up_6_3_opt + subroutine coeffs_up_6 ( a1, q1, a2, q2 ) + + CCTK_REAL, dimension(-4:4), intent(OUT) :: a1, a2 + CCTK_REAL, dimension(5,4), intent(OUT) :: q1, q2 + + a1(-4) = 0.01666666666666666666666666666666666666667_wp + a1(-3) = -0.1333333333333333333333333333333333333333_wp + a1(-2) = 0.5_wp + a1(-1) = -1.333333333333333333333333333333333333333_wp + a1(0) = 0.5833333333333333333333333333333333333333_wp + a1(1) = 0.4_wp + a1(2) = -0.03333333333333333333333333333333333333333_wp + a1(3) = zero + a1(4) = zero + + q1(1,1) = zero + q1(2,1) = zero + q1(3,1) = zero + q1(4,1) = zero + q1(5,1) = zero + + q1(1,2) = -1.0_wp + q1(2,2) = 1.0_wp + q1(3,2) = zero + q1(4,2) = zero + q1(5,2) = zero + + q1(1,3) = 0.5_wp + q1(2,3) = -2.0_wp + q1(3,3) = 1.5_wp + q1(4,3) = zero + q1(5,3) = zero + + q1(1,4) = -0.08333333333333333333333333333333333333333_wp + q1(2,4) = 0.5_wp + q1(3,4) = -1.5_wp + q1(4,4) = 0.8333333333333333333333333333333333333333_wp + q1(5,4) = 0.25_wp + + a2(-4) = zero + a2(-3) = zero + a2(-2) = 0.03333333333333333333333333333333333333333_wp + a2(-1) = -0.4_wp + a2(0) = -0.5833333333333333333333333333333333333333_wp + a2(1) = 1.333333333333333333333333333333333333333_wp + a2(2) = -0.5_wp + a2(3) = 0.1333333333333333333333333333333333333333_wp + a2(4) = -0.01666666666666666666666666666666666666667_wp + + q2(1,1) = zero + q2(2,1) = zero + q2(3,1) = zero + q2(4,1) = zero + q2(5,1) = zero + + q2(1,2) = zero + q2(2,2) = zero + q2(3,2) = zero + q2(4,2) = -1.0_wp + q2(5,2) = 1.0_wp + + q2(1,3) = zero + q2(2,3) = zero + q2(3,3) = -1.5_wp + q2(4,3) = 2.0_wp + q2(5,3) = -0.5_wp + + q2(1,4) = -0.25_wp + q2(2,4) = -0.8333333333333333333333333333333333333333_wp + q2(3,4) = 1.5_wp + q2(4,4) = -0.5_wp + q2(5,4) = 0.08333333333333333333333333333333333333333_wp + + end subroutine coeffs_up_6 + subroutine coeffs_2_6_3 ( a, q ) CCTK_REAL, dimension(4), intent(OUT) :: a @@ -1201,6 +1364,116 @@ module All_Coeffs_mod end subroutine coeffs_up_8_4_opt + subroutine coeffs_up_8 ( a1, q1, a2, q2 ) + + CCTK_REAL, dimension(-5:5), intent(OUT) :: a1, a2 + CCTK_REAL, dimension(7,5), intent(OUT) :: q1, q2 + + a1(-5) = -0.003571428571428571428571428571428571428571_wp + a1(-4) = 0.03571428571428571428571428571428571428571_wp + a1(-3) = -0.1666666666666666666666666666666666666667_wp + a1(-2) = 0.5_wp + a1(-1) = -1.25_wp + a1(0) = 0.45_wp + a1(1) = 0.5_wp + a1(2) = -0.07142857142857142857142857142857142857143_wp + a1(3) = 0.005952380952380952380952380952380952380952_wp + a1(4) = zero + a1(5) = zero + + q1(2,1) = zero + q1(3,1) = zero + q1(4,1) = zero + q1(5,1) = zero + q1(6,1) = zero + q1(7,1) = zero + + q1(1,2) = -1.0_wp + q1(2,2) = 1.0_wp + q1(3,2) = zero + q1(4,2) = zero + q1(5,2) = zero + q1(6,2) = zero + q1(7,2) = zero + + q1(1,3) = 0.5_wp + q1(2,3) = -2.0_wp + q1(3,3) = 1.5_wp + q1(4,3) = zero + q1(5,3) = zero + q1(6,3) = zero + q1(7,3) = zero + + q1(1,4) = -0.08333333333333333333333333333333333333333_wp + q1(2,4) = 0.5_wp + q1(3,4) = -1.5_wp + q1(4,4) = 0.8333333333333333333333333333333333333333_wp + q1(5,4) = 0.25_wp + q1(6,4) = zero + q1(7,4) = zero + + q1(1,5) = 0.01666666666666666666666666666666666666667_wp + q1(2,5) = -0.1333333333333333333333333333333333333333_wp + q1(3,5) = 0.5_wp + q1(4,5) = -1.333333333333333333333333333333333333333_wp + q1(5,5) = 0.5833333333333333333333333333333333333333_wp + q1(6,5) = 0.4_wp + q1(7,5) = -0.03333333333333333333333333333333333333333_wp + + a2(-5) = zero + a2(-4) = zero + a2(-3) = -0.005952380952380952380952380952380952380952_wp + a2(-2) = 0.07142857142857142857142857142857142857143_wp + a2(-1) = -0.5_wp + a2(0) = -0.45_wp + a2(1) = 1.25_wp + a2(2) = -0.5_wp + a2(3) = 0.1666666666666666666666666666666666666667_wp + a2(4) = -0.03571428571428571428571428571428571428571_wp + a2(5) = 0.003571428571428571428571428571428571428571_wp + + q2(1,1) = zero + q2(2,1) = zero + q2(3,1) = zero + q2(4,1) = zero + q2(5,1) = zero + q2(6,1) = zero + q2(7,1) = zero + + q2(1,2) = zero + q2(2,2) = zero + q2(3,2) = zero + q2(4,2) = zero + q2(5,2) = zero + q2(6,2) = -1.0_wp + q2(7,2) = 1.0_wp + + q2(1,3) = zero + q2(2,3) = zero + q2(3,3) = zero + q2(4,3) = zero + q2(5,3) = -1.5_wp + q2(6,3) = 2.0_wp + q2(7,3) = -0.5_wp + + q2(1,4) = zero + q2(2,4) = zero + q2(3,4) = -0.25_wp + q2(4,4) = -0.8333333333333333333333333333333333333333_wp + q2(5,4) = 1.5_wp + q2(6,4) = -0.5_wp + q2(7,4) = 0.08333333333333333333333333333333333333333_wp + + q2(1,5) = 0.03333333333333333333333333333333333333333_wp + q2(2,5) = -0.4_wp + q2(3,5) = -0.5833333333333333333333333333333333333333_wp + q2(4,5) = 1.333333333333333333333333333333333333333_wp + q2(5,5) = -0.5_wp + q2(6,5) = 0.1333333333333333333333333333333333333333_wp + q2(7,5) = -0.01666666666666666666666666666666666666667_wp + + end subroutine coeffs_up_8 + subroutine coeffs_2_8_4 ( a, q ) CCTK_REAL, dimension(5), intent(OUT) :: a diff --git a/src/Coefficients_2_1.F90 b/src/Coefficients_2_1.F90 index 99f21c7..3f89b83 100644 --- a/src/Coefficients_2_1.F90 +++ b/src/Coefficients_2_1.F90 @@ -136,3 +136,85 @@ subroutine set_coeff_up_2_1 ( nsize, dir, loc_order, bb, gsize, imin, imax, dd ) end if end do end subroutine set_coeff_up_2_1 + +subroutine set_coeff_up_2 ( nsize, dir, loc_order, bb, gsize, imin, imax, dd ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(IN) :: nsize, dir, loc_order + CCTK_INT, dimension(2), intent(IN) :: bb + CCTK_INT, intent(IN) :: gsize + CCTK_INT, dimension(nsize), intent(OUT) :: imin, imax + CCTK_REAL, dimension(nsize,nsize), intent(OUT) :: dd + + CCTK_REAL, dimension(-2:2), save :: a1, a2 + CCTK_REAL, dimension(3,2), save :: q1, q2 + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( .not. ( dir == -1 .or. dir == 1 ) ) then + call CCTK_WARN(0, 'Internal error: set_coeff_up_2 called with invalid dir value') + end if + + if ( first ) then + call coeffs_up_2 ( a1, q1, a2, q2 ) + first = .false. + end if + + dd = zero + imin = 0 + imax = -1 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + if ( dir == -1 ) then + dd(1,1) = q1(1,1); + imin(1) = 1; imax(1) = 1 + + dd(1:2,2) = q1(1:2,2); + imin(2) = 1; imax(2) = 2 + else + dd(1,1) = q2(1,1); + imin(1) = 1; imax(1) = 1 + + dd(2:4,2) = a2(0:2); + imin(2) = 2; imax(2) = 4 + end if + il = 3 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + if ( dir == -1 ) then + dd(nsize,nsize) = q1(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + dd(nsize-3:nsize-1,nsize-1) = a1(-2:0) + imin(nsize-1) = nsize-3; imax(nsize-1) = nsize-1 + else + dd(nsize,nsize) = q2(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + dd(nsize-1:nsize,nsize-1) = q2(2:3,2) + imin(nsize-1) = nsize-1; imax(nsize-1) = nsize + end if + ir = nsize - 2 + end if + do i = il, ir + if ( dir == -1 ) then + dd(i-2:i,i) = a1(-2:0) + imin(i) = i-2; imax(i) = i + else + dd(i:i+2,i) = a2(0:2) + imin(i) = i; imax(i) = i+2 + end if + end do +end subroutine set_coeff_up_2 diff --git a/src/Coefficients_4_2.F90 b/src/Coefficients_4_2.F90 index 099c5b6..d2b1b07 100644 --- a/src/Coefficients_4_2.F90 +++ b/src/Coefficients_4_2.F90 @@ -178,3 +178,97 @@ subroutine set_coeff_up_4_2 ( nsize, dir, loc_order, bb, gsize, imin, imax, dd ) end if end do end subroutine set_coeff_up_4_2 + +subroutine set_coeff_up_4 ( nsize, dir, loc_order, bb, gsize, imin, imax, dd ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(IN) :: nsize, dir, loc_order + CCTK_INT, dimension(2), intent(IN) :: bb + CCTK_INT, intent(IN) :: gsize + CCTK_INT, dimension(nsize), intent(OUT) :: imin, imax + CCTK_REAL, dimension(nsize,nsize), intent(OUT) :: dd + + CCTK_REAL, dimension(-3:3), save :: a1, a2 + CCTK_REAL, dimension(4,3), save :: q1, q2 + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( .not. ( dir == -1 .or. dir == 1 ) ) then + call CCTK_WARN(0, 'Internal error: set_coeff_up_2 called with invalid dir value') + end if + + if ( first ) then + call coeffs_up_4 ( a1, q1, a2, q2 ) + first = .false. + end if + + dd = zero + imin = 0 + imax = -1 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + if ( dir == -1 ) then + dd(1,1) = q1(1,1); + imin(1) = 1; imax(1) = 1 + + dd(1:2,2) = q1(1:2,2); + imin(2) = 1; imax(2) = 2 + + dd(1:3,3) = q1(1:3,3); + imin(3) = 1; imax(3) = 3 + else + dd(1,1) = q2(1,1); + imin(1) = 1; imax(1) = 1 + + dd(1:5,2) = a2(-1:3); + imin(2) = 1; imax(2) = 5 + + dd(2:6,3) = a2(-1:3); + imin(3) = 2; imax(3) = 6 + end if + il = 4 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + if ( dir == -1 ) then + dd(nsize,nsize) = q1(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + dd(nsize-4:nsize,nsize-1) = a1(-3:1) + imin(nsize-1) = nsize-4; imax(nsize-1) = nsize + + dd(nsize-5:nsize-1,nsize-2) = a1(-3:1) + imin(nsize-2) = nsize-5; imax(nsize-2) = nsize-1 + else + dd(nsize,nsize) = q2(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + dd(nsize-1:nsize,nsize-1) = q2(3:4,2) + imin(nsize-1) = nsize-1; imax(nsize-1) = nsize + + dd(nsize-2:nsize,nsize-2) = q2(2:4,3) + imin(nsize-2) = nsize-2; imax(nsize-2) = nsize + end if + ir = nsize - 3 + end if + do i = il, ir + if ( dir == -1 ) then + dd(i-3:i+1,i) = a1(-3:1) + imin(i) = i-3; imax(i) = i+1 + else + dd(i-1:i+3,i) = a2(-1:3) + imin(i) = i-1; imax(i) = i+3 + end if + end do +end subroutine set_coeff_up_4 diff --git a/src/Coefficients_6_3_min_err_coeff.F90 b/src/Coefficients_6_3_min_err_coeff.F90 index c0eef1e..04d4e4f 100644 --- a/src/Coefficients_6_3_min_err_coeff.F90 +++ b/src/Coefficients_6_3_min_err_coeff.F90 @@ -215,3 +215,109 @@ subroutine set_coeff_up_6_3_opt ( nsize, dir, loc_order, bb, gsize, imin, imax, end if end do end subroutine set_coeff_up_6_3_opt + +subroutine set_coeff_up_6 ( nsize, dir, loc_order, bb, gsize, imin, imax, dd ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(IN) :: nsize, dir, loc_order + CCTK_INT, dimension(2), intent(IN) :: bb + CCTK_INT, intent(IN) :: gsize + CCTK_INT, dimension(nsize), intent(OUT) :: imin, imax + CCTK_REAL, dimension(nsize,nsize), intent(OUT) :: dd + + CCTK_REAL, dimension(-4:4), save :: a1, a2 + CCTK_REAL, dimension(5,4), save :: q1, q2 + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( .not. ( dir == -1 .or. dir == 1 ) ) then + call CCTK_WARN(0, 'Internal error: set_coeff_up_2 called with invalid dir value') + end if + + if ( first ) then + call coeffs_up_6 ( a1, q1, a2, q2 ) + first = .false. + end if + + dd = zero + imin = 0 + imax = -1 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + if ( dir == -1 ) then + dd(1,1) = q1(1,1); + imin(1) = 1; imax(1) = 1 + + dd(1:2,2) = q1(1:2,2); + imin(2) = 1; imax(2) = 2 + + dd(1:3,3) = q1(1:3,3); + imin(3) = 1; imax(3) = 3 + + dd(1:5,4) = q1(1:5,4); + imin(4) = 1; imax(4) = 5 + else + dd(1,1) = q2(1,1); + imin(1) = 1; imax(1) = 1 + + dd(1:5,2) = q2(1:5,4); + imin(2) = 1; imax(2) = 5 + + dd(1:7,3) = a2(-2:4); + imin(3) = 1; imax(3) = 7 + + dd(2:8,4) = a2(-2:4); + imin(4) = 2; imax(4) = 8 + end if + il = 5 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + if ( dir == -1 ) then + dd(nsize,nsize) = q1(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + dd(nsize-4:nsize,nsize-1) = q1(1:5,4) + imin(nsize-1) = nsize-4; imax(nsize-1) = nsize + + dd(nsize-6:nsize,nsize-2) = a1(-4:2) + imin(nsize-2) = nsize-6; imax(nsize-2) = nsize + + dd(nsize-7:nsize-1,nsize-3) = a1(-4:2) + imin(nsize-3) = nsize-7; imax(nsize-3) = nsize-1 + else + dd(nsize,nsize) = q2(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + dd(nsize-1:nsize,nsize-1) = q2(4:5,2) + imin(nsize-1) = nsize-1; imax(nsize-1) = nsize + + dd(nsize-2:nsize,nsize-2) = q2(3:5,3) + imin(nsize-2) = nsize-2; imax(nsize-2) = nsize + + dd(nsize-4:nsize,nsize-3) = q2(1:5,4) + imin(nsize-3) = nsize-4; imax(nsize-3) = nsize + end if + ir = nsize - 4 + end if + do i = il, ir + if ( dir == -1 ) then + dd(i-4:i+2,i) = a1(-4:2) + imin(i) = i-4; imax(i) = i+2 + else + dd(i-2:i+4,i) = a2(-2:4) + imin(i) = i-2; imax(i) = i+4 + end if + end do +end subroutine set_coeff_up_6 diff --git a/src/Coefficients_8_4_min_err_coeff.F90 b/src/Coefficients_8_4_min_err_coeff.F90 index 5bae502..8f0c5bf 100644 --- a/src/Coefficients_8_4_min_err_coeff.F90 +++ b/src/Coefficients_8_4_min_err_coeff.F90 @@ -251,3 +251,121 @@ subroutine set_coeff_up_8_4_opt ( nsize, dir, loc_order, bb, gsize, imin, imax, end if end do end subroutine set_coeff_up_8_4_opt + +subroutine set_coeff_up_8 ( nsize, dir, loc_order, bb, gsize, imin, imax, dd ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(IN) :: nsize, dir, loc_order + CCTK_INT, dimension(2), intent(IN) :: bb + CCTK_INT, intent(IN) :: gsize + CCTK_INT, dimension(nsize), intent(OUT) :: imin, imax + CCTK_REAL, dimension(nsize,nsize), intent(OUT) :: dd + + CCTK_REAL, dimension(-5:5), save :: a1, a2 + CCTK_REAL, dimension(7,5), save :: q1, q2 + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( .not. ( dir == -1 .or. dir == 1 ) ) then + call CCTK_WARN(0, 'Internal error: set_coeff_up_8 called with invalid dir value') + end if + + if ( first ) then + call coeffs_up_8 ( a1, q1, a2, q2 ) + first = .false. + end if + + dd = zero + imin = 0 + imax = -1 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + if ( dir == -1 ) then + dd(1,1) = q1(1,1); + imin(1) = 1; imax(1) = 1 + + dd(1:2,2) = q1(1:2,2); + imin(2) = 1; imax(2) = 2 + + dd(1:3,3) = q1(1:3,3); + imin(3) = 1; imax(3) = 3 + + dd(1:5,4) = q1(1:5,4); + imin(4) = 1; imax(4) = 5 + + dd(1:7,5) = q1(1:7,5); + imin(5) = 1; imax(5) = 7 + else + dd(1,1) = q2(1,1); + imin(1) = 1; imax(1) = 1 + + dd(1:5,2) = q2(3:7,4); + imin(2) = 1; imax(2) = 5 + + dd(1:7,3) = q2(1:7,5); + imin(3) = 1; imax(3) = 7 + + dd(1:9,4) = a2(-3:5); + imin(4) = 1; imax(4) = 9 + + dd(2:10,5) = a2(-3:5); + imin(5) = 2; imax(5) = 10 + end if + il = 6 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + if ( dir == -1 ) then + dd(nsize,nsize) = q1(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + dd(nsize-4:nsize,nsize-1) = q1(1:5,4) + imin(nsize-1) = nsize-4; imax(nsize-1) = nsize + + dd(nsize-6:nsize,nsize-2) = q1(1:7,5) + imin(nsize-2) = nsize-6; imax(nsize-2) = nsize + + dd(nsize-8:nsize,nsize-3) = a1(-5:3) + imin(nsize-3) = nsize-8; imax(nsize-3) = nsize + + dd(nsize-9:nsize-1,nsize-4) = a1(-5:3) + imin(nsize-4) = nsize-9; imax(nsize-4) = nsize-1 + else + dd(nsize,nsize) = q2(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + dd(nsize-1:nsize,nsize-1) = q2(6:7,2) + imin(nsize-1) = nsize-1; imax(nsize-1) = nsize + + dd(nsize-2:nsize,nsize-2) = q2(5:7,3) + imin(nsize-2) = nsize-2; imax(nsize-2) = nsize + + dd(nsize-4:nsize,nsize-3) = q2(3:7,4) + imin(nsize-3) = nsize-4; imax(nsize-3) = nsize + + dd(nsize-6:nsize,nsize-4) = q2(1:7,5) + imin(nsize-4) = nsize-6; imax(nsize-4) = nsize + end if + ir = nsize - 5 + end if + do i = il, ir + if ( dir == -1 ) then + dd(i-5:i+3,i) = a1(-5:3) + imin(i) = i-5; imax(i) = i+3 + else + dd(i-3:i+5,i) = a2(-3:5) + imin(i) = i-3; imax(i) = i+5 + end if + end do +end subroutine set_coeff_up_8 diff --git a/src/get_up_coeffs.c b/src/get_up_coeffs.c index 3a78f6a..e75ab9a 100644 --- a/src/get_up_coeffs.c +++ b/src/get_up_coeffs.c @@ -60,6 +60,39 @@ void DiffUpCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_INT *imax, CCTK_REAL *q); + void CCTK_FCALL CCTK_FNAME(set_coeff_up_2)(const CCTK_INT *nsize, + const CCTK_INT *up_dir, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); + void CCTK_FCALL CCTK_FNAME(set_coeff_up_4)(const CCTK_INT *nsize, + const CCTK_INT *up_dir, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); + void CCTK_FCALL CCTK_FNAME(set_coeff_up_6)(const CCTK_INT *nsize, + const CCTK_INT *up_dir, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); + void CCTK_FCALL CCTK_FNAME(set_coeff_up_8)(const CCTK_INT *nsize, + const CCTK_INT *up_dir, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); + ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2]; if ( table_handle >=0 ) { @@ -99,35 +132,58 @@ void DiffUpCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, } if ( CCTK_Equals(norm_type,"Diagonal") ) { - if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) { - switch(loc_order) { - case 2: { - CCTK_FNAME(set_coeff_up_2_1)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); - break; - } - case 4: { - CCTK_FNAME(set_coeff_up_4_2)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); - break; - } - default: - CCTK_WARN (0, "Unknown upwinding stencil specified"); + if ( sbp_upwind_deriv ) { + if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) { + switch(loc_order) { + case 2: { + CCTK_FNAME(set_coeff_up_2_1)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 4: { + CCTK_FNAME(set_coeff_up_4_2)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + break; + } + default: + CCTK_WARN (0, "Unknown upwinding stencil specified"); + } + } else { + switch(loc_order) { + case 2: { + CCTK_FNAME(set_coeff_up_2_1)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 4: { + CCTK_FNAME(set_coeff_up_4_2)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 6: { + CCTK_FNAME(set_coeff_up_6_3_opt)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 8: { + CCTK_FNAME(set_coeff_up_8_4_opt)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + break; + } + default: + CCTK_WARN (0, "Unknown upwinding stencil specified"); + } } - } else { + } else { switch(loc_order) { case 2: { - CCTK_FNAME(set_coeff_up_2_1)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff_up_2)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); break; } case 4: { - CCTK_FNAME(set_coeff_up_4_2)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff_up_4)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); break; } case 6: { - CCTK_FNAME(set_coeff_up_6_3_opt)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff_up_6)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); break; } case 8: { - CCTK_FNAME(set_coeff_up_8_4_opt)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff_up_8)(&nsize,&up,&loc_order,bb,&gsize,imin,imax,q); break; } default: -- cgit v1.2.3