From 9fbf1f753b7927be6d888b46049953f36a1fe3f4 Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 31 Mar 2009 23:00:16 +0000 Subject: Introduce a boolean parameter to switch of SBP near outer boundary. Instead switch to progressively lower order centered stencil as the boundery is approached, similar to what is done with the second order and upwind stencils. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@118 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- param.ccl | 4 ++ src/All_Coeffs_mod.F90 | 102 ++++++++++++++++++++++++++++++++++++++++++++++++ src/Coefficients_2.F90 | 56 +++++++++++++++++++++++++++ src/Coefficients_4.F90 | 63 ++++++++++++++++++++++++++++++ src/Coefficients_6.F90 | 69 +++++++++++++++++++++++++++++++++ src/Coefficients_8.F90 | 75 +++++++++++++++++++++++++++++++++++ src/get_coeffs.c | 103 ++++++++++++++++++++++++++++++++++++------------- src/make.code.defn | 4 ++ 8 files changed, 450 insertions(+), 26 deletions(-) create mode 100644 src/Coefficients_2.F90 create mode 100644 src/Coefficients_4.F90 create mode 100644 src/Coefficients_6.F90 create mode 100644 src/Coefficients_8.F90 diff --git a/param.ccl b/param.ccl index d4314cf..d382d36 100644 --- a/param.ccl +++ b/param.ccl @@ -15,6 +15,10 @@ KEYWORD operator_type "Type of operator" STEERABLE=always "Optimized" :: "Optimized for performance" } "Optimized" +BOOLEAN sbp_1st_deriv "Should the 1st derivative operator be SBP" STEERABLE=always +{ +} "yes" + BOOLEAN sbp_2nd_deriv "Should the 2nd derivative operator be SBP" STEERABLE=always { } "yes" diff --git a/src/All_Coeffs_mod.F90 b/src/All_Coeffs_mod.F90 index 1886497..55a76ab 100644 --- a/src/All_Coeffs_mod.F90 +++ b/src/All_Coeffs_mod.F90 @@ -28,6 +28,17 @@ module All_Coeffs_mod end subroutine coeffs_1_2_1 + subroutine coeffs_1_2 ( a, q ) + + CCTK_REAL, dimension(1), intent(OUT) :: a + CCTK_REAL, dimension(1,1), intent(OUT) :: q + + a(1) = 1.0_wp/2.0_wp + + q(1,1) = zero + + end subroutine coeffs_1_2 + subroutine coeffs_up_2_1 ( a1, q1, a2, q2 ) CCTK_REAL, dimension(-1:1), intent(OUT) :: a1, a2 @@ -161,6 +172,24 @@ module All_Coeffs_mod end subroutine coeffs_1_4_2 + subroutine coeffs_1_4 ( a, q ) + + CCTK_REAL, dimension(2), intent(OUT) :: a + CCTK_REAL, dimension(3,2), intent(OUT) :: q + + a(1) = 2.0_wp/3.0_wp + a(2) = -1.0_wp/12.0_wp + + q(1,1) = zero + q(2,1) = zero + q(3,1) = zero + + q(1,2) = -1.0_wp/2.0_wp + q(2,2) = zero + q(3,2) = 1.0_wp/2.0_wp + + end subroutine coeffs_1_4 + subroutine coeffs_up_4_2 ( a1, q1, a2, q2 ) CCTK_REAL, dimension(-2:2), intent(OUT) :: a1, a2 @@ -526,6 +555,35 @@ module All_Coeffs_mod end subroutine coeffs_1_6_3_opt + subroutine coeffs_1_6 ( a, q ) + + CCTK_REAL, dimension(3), intent(OUT) :: a + CCTK_REAL, dimension(5,3), intent(OUT) :: q + + a(1) = 3.0_wp/4.0_wp + a(2) = -3.0_wp/20.0_wp + a(3) = 1.0_wp/60.0_wp + + q(1,1) = zero + q(2,1) = zero + q(3,1) = zero + q(4,1) = zero + q(5,1) = zero + + q(1,2) = -1.0_wp/2.0_wp + q(2,2) = zero + q(3,2) = 1.0_wp/2.0_wp + q(4,2) = zero + q(5,2) = zero + + q(1,3) = 1.0_wp/12.0_wp + q(2,3) = -2.0_wp/3.0_wp + q(3,3) = zero + q(4,3) = 2.0_wp/3.0_wp + q(5,3) = -1.0_wp/12.0_wp + + end subroutine coeffs_1_6 + subroutine coeffs_up_6_3_opt ( a1, q1, a2, q2 ) CCTK_REAL, dimension(-3:3), intent(OUT) :: a1, a2 @@ -1129,6 +1187,50 @@ module All_Coeffs_mod end subroutine coeffs_1_8_4_opt + subroutine coeffs_1_8 ( a, q ) + + CCTK_REAL, dimension(4), intent(OUT) :: a + CCTK_REAL, dimension(7,4), intent(OUT) :: q + + a(1) = 4.0_wp/5.0_wp + a(2) = -1.0_wp/5.0_wp + a(3) = 4.0_wp/105.0_wp + a(4) = -1.0_wp/280.0_wp + + q(1,1) = zero + q(2,1) = zero + q(3,1) = zero + q(4,1) = zero + q(5,1) = zero + q(6,1) = zero + q(7,1) = zero + + q(1,2) = -1.0_wp/2.0_wp + q(2,2) = zero + q(3,2) = 1.0_wp/2.0_wp + q(4,2) = zero + q(5,2) = zero + q(6,2) = zero + q(7,1) = zero + + q(1,3) = 1.0_wp/12.0_wp + q(2,3) = -2.0_wp/3.0_wp + q(3,3) = zero + q(4,3) = 2.0_wp/3.0_wp + q(5,3) = -1.0_wp/12.0_wp + q(6,3) = zero + q(7,3) = zero + + q(1,4) = -1.0_wp/60.0_wp + q(2,4) = 3.0_wp/20.0_wp + q(3,4) = -3.0_wp/4.0_wp + q(4,4) = zero + q(5,4) = 3.0_wp/4.0_wp + q(6,4) = -3.0_wp/20.0_wp + q(7,4) = 1.0_wp/60.0_wp + + end subroutine coeffs_1_8 + subroutine coeffs_up_8_4_opt ( a1, q1, a2, q2 ) CCTK_REAL, dimension(-4:4), intent(OUT) :: a1, a2 diff --git a/src/Coefficients_2.F90 b/src/Coefficients_2.F90 new file mode 100644 index 0000000..122961b --- /dev/null +++ b/src/Coefficients_2.F90 @@ -0,0 +1,56 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine set_coeff_2 ( nsize, loc_order, bb, gsize, imin, imax, dd ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(IN) :: nsize, 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(1), save :: a + CCTK_REAL, dimension(1,1), save :: q + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( first ) then + call coeffs_1_2 ( a, q ) + first = .false. + end if + + dd = zero + imin = 0 + imax = -1 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dd(1,1) = q(1,1); + imin(1) = 1; imax(1) = 1 + il = 2 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + dd(nsize,nsize) = q(1,1) + imin(nsize) = nsize; imax(nsize) = nsize + ir = nsize - 1 + end if + do i = il, ir + dd(i-1,i) = -a(1); + dd(i+1,i) = a(1) + imin(i) = i-1; imax(i) = i+1 + end do + +end subroutine set_coeff_2 diff --git a/src/Coefficients_4.F90 b/src/Coefficients_4.F90 new file mode 100644 index 0000000..2211559 --- /dev/null +++ b/src/Coefficients_4.F90 @@ -0,0 +1,63 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine set_coeff_4 ( nsize, loc_order, bb, gsize, imin, imax, dd ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(IN) :: nsize, 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), save :: a + CCTK_REAL, dimension(3,2), save :: q + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( first ) then + call coeffs_1_4 ( a, q ) + first = .false. + end if + + dd = zero + imin = 0 + imax = -1 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dd(1:3,1) = q(1:3,1) + imin(1) = 1; imax(1) = 1 + + dd(1:3,2) = q(1:3,2) + imin(2) = 1; imax(2) = 3 + + il = 3 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + dd(nsize-2:nsize,nsize-1) = -q(3:1:-1,2) + imin(nsize-1) = nsize-2; imax(nsize-1) = nsize + + dd(nsize-2:nsize,nsize) = -q(3:1:-1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + ir = nsize - 2 + end if + do i = il, ir + dd(i-2:i-1,i) = -a(2:1:-1); + dd(i+1:i+2,i) = a(1:2) + imin(i) = i-2; imax(i) = i+2 + end do +end subroutine set_coeff_4 diff --git a/src/Coefficients_6.F90 b/src/Coefficients_6.F90 new file mode 100644 index 0000000..a27089f --- /dev/null +++ b/src/Coefficients_6.F90 @@ -0,0 +1,69 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine set_coeff_6 ( nsize, loc_order, bb, gsize, imin, imax, dd ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(IN) :: nsize, 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), save :: a + CCTK_REAL, dimension(5,3), save :: q + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( first ) then + call coeffs_1_6 ( a, q ) + first = .false. + end if + + dd = zero + imin = 0 + imax = -1 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dd(1:5,1) = q(1:5,1) + imin(1) = 1; imax(1) = 1 + + dd(1:5,2) = q(1:5,2) + imin(2) = 1; imax(2) = 3 + + dd(1:5,3) = q(1:5,3) + imin(3) = 1; imax(3) = 5 + + il = 4 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + dd(nsize-4:nsize,nsize-2) = -q(5:1:-1,3) + imin(nsize-2) = nsize-4; imax(nsize-2) = nsize + + dd(nsize-4:nsize,nsize-1) = -q(5:1:-1,2) + imin(nsize-1) = nsize-2; imax(nsize-1) = nsize + + dd(nsize-4:nsize,nsize) = -q(5:1:-1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + ir = nsize - 3 + end if + do i = il, ir + dd(i-3:i-1,i) = -a(3:1:-1); + dd(i+1:i+3,i) = a(1:3) + imin(i) = i-3; imax(i) = i+3 + end do +end subroutine set_coeff_6 diff --git a/src/Coefficients_8.F90 b/src/Coefficients_8.F90 new file mode 100644 index 0000000..a9fdeb6 --- /dev/null +++ b/src/Coefficients_8.F90 @@ -0,0 +1,75 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine set_coeff_8 ( nsize, loc_order, bb, gsize, imin, imax, dd ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(IN) :: nsize, 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), save :: a + CCTK_REAL, dimension(7,4), save :: q + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( first ) then + call coeffs_1_8 ( a, q ) + first = .false. + end if + + dd = zero + imin = 0 + imax = -1 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dd(1:7,1) = q(1:7,1) + imin(1) = 1; imax(1) = 1 + + dd(1:7,2) = q(1:7,2) + imin(2) = 1; imax(2) = 3 + + dd(1:7,3) = q(1:7,3) + imin(3) = 1; imax(3) = 5 + + dd(1:7,4) = q(1:7,4) + imin(4) = 1; imax(4) = 7 + + il = 5 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + dd(nsize-6:nsize,nsize-3) = -q(7:1:-1,4) + imin(nsize-3) = nsize-6; imax(nsize-3) = nsize + + dd(nsize-6:nsize,nsize-2) = -q(7:1:-1,3) + imin(nsize-2) = nsize-4; imax(nsize-2) = nsize + + dd(nsize-6:nsize,nsize-1) = -q(7:1:-1,2) + imin(nsize-1) = nsize-2; imax(nsize-1) = nsize + + dd(nsize-6:nsize,nsize) = -q(7:1:-1,1) + imin(nsize) = nsize; imax(nsize) = nsize + + ir = nsize - 4 + end if + do i = il, ir + dd(i-4:i-1,i) = -a(4:1:-1); + dd(i+1:i+4,i) = a(1:4) + imin(i) = i-4; imax(i) = i+4 + end do +end subroutine set_coeff_8 diff --git a/src/get_coeffs.c b/src/get_coeffs.c index 836ff52..31b3601 100644 --- a/src/get_coeffs.c +++ b/src/get_coeffs.c @@ -32,6 +32,13 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_INT *imin, CCTK_INT *imax, CCTK_REAL *q); + void CCTK_FCALL CCTK_FNAME(set_coeff_2)(const CCTK_INT *nsize, + 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_4_2)(const CCTK_INT *nsize, const CCTK_INT *loc_order, const CCTK_INT *bb, @@ -39,6 +46,13 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_INT *imin, CCTK_INT *imax, CCTK_REAL *q); + void CCTK_FCALL CCTK_FNAME(set_coeff_4)(const CCTK_INT *nsize, + 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_4_3)(const CCTK_INT *nsize, const CCTK_INT *loc_order, const CCTK_INT *bb, @@ -67,6 +81,13 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_INT *imin, CCTK_INT *imax, CCTK_REAL *q); + void CCTK_FCALL CCTK_FNAME(set_coeff_6)(const CCTK_INT *nsize, + 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_6_5)(const CCTK_INT *nsize, const CCTK_INT *loc_order, const CCTK_INT *bb, @@ -95,6 +116,13 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_INT *imin, CCTK_INT *imax, CCTK_REAL *q); + void CCTK_FCALL CCTK_FNAME(set_coeff_8)(const CCTK_INT *nsize, + 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]; @@ -135,47 +163,70 @@ void DiffCoeff ( 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_2_1)(&nsize,&loc_order,bb,&gsize,imin,imax,q); - break; - } - case 4: { - CCTK_FNAME(set_coeff_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q); - break; - } - case 6: { - CCTK_FNAME(set_coeff_6_3)(&nsize,&loc_order,bb,&gsize,imin,imax,q); - break; - } - case 8: { - CCTK_FNAME(set_coeff_8_4)(&nsize,&loc_order,bb,&gsize,imin,imax,q); - break; - } - default: - CCTK_WARN (0, "Unknown stencil specified"); + if ( sbp_1st_deriv ) { + if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) { + switch(loc_order) { + case 2: { + CCTK_FNAME(set_coeff_2_1)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 4: { + CCTK_FNAME(set_coeff_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 6: { + CCTK_FNAME(set_coeff_6_3)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 8: { + CCTK_FNAME(set_coeff_8_4)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } + default: + CCTK_WARN (0, "Unknown stencil specified"); + } + } else { + switch(loc_order) { + case 2: { + CCTK_FNAME(set_coeff_2_1)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 4: { + CCTK_FNAME(set_coeff_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 6: { + CCTK_FNAME(set_coeff_6_3_opt)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } + case 8: { + CCTK_FNAME(set_coeff_8_4_opt)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } + default: + CCTK_WARN (0, "Unknown stencil specified"); + } } } else { - switch(loc_order) { + switch(loc_order) { case 2: { - CCTK_FNAME(set_coeff_2_1)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q); break; } case 4: { - CCTK_FNAME(set_coeff_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff_4)(&nsize,&loc_order,bb,&gsize,imin,imax,q); break; } case 6: { - CCTK_FNAME(set_coeff_6_3_opt)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff_6)(&nsize,&loc_order,bb,&gsize,imin,imax,q); break; } case 8: { - CCTK_FNAME(set_coeff_8_4_opt)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff_8)(&nsize,&loc_order,bb,&gsize,imin,imax,q); break; } default: - CCTK_WARN (0, "Unknown stencil specified"); + CCTK_WARN (0, "Unknown 1st derivative stencil specified"); } } } else { diff --git a/src/make.code.defn b/src/make.code.defn index 8eaaf16..7fa5829 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -46,19 +46,23 @@ SRCS = call_derivs.c \ get_up_coeffs.c \ get_coeffs2.c \ Coefficients_2_1.F90 \ + Coefficients_2.F90 \ Coefficients2_2_1.F90 \ Coefficients2_2.F90 \ Coefficients_4_2.F90 \ + Coefficients_4.F90 \ Coefficients2_4_2.F90 \ Coefficients2_4_2_min_err_coeff.F90 \ Coefficients2_4.F90 \ Coefficients_6_3.F90 \ Coefficients2_6_3.F90 \ Coefficients_6_3_min_err_coeff.F90 \ + Coefficients_6.F90 \ Coefficients2_6.F90 \ Coefficients_8_4.F90 \ Coefficients2_8_4.F90 \ Coefficients_8_4_min_err_coeff.F90 \ + Coefficients_8.F90 \ Coefficients2_8.F90 \ Coefficients_4_3.F90 \ Coefficients_4_3_min_err_coeff.F90 \ -- cgit v1.2.3