From 2a45c231d2a5b7a1e29faaf599f3a19cf91b9293 Mon Sep 17 00:00:00 2001 From: diener Date: Thu, 18 Oct 2007 08:58:37 +0000 Subject: The second derivative in the 4-2 case is not unique, so I added the minimal bandwidth operator in addition to the ABTE operator already there. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@86 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- src/All_Coeffs_mod.F90 | 41 +++++++- src/Derivatives2_4_2.F90 | 8 +- src/Derivatives2_4_2_min_err_coeff.F90 | 177 +++++++++++++++++++++++++++++++++ src/call_derivs2.c | 44 ++++++-- src/make.code.defn | 1 + src/make.code.deps | 1 + 6 files changed, 257 insertions(+), 15 deletions(-) create mode 100644 src/Derivatives2_4_2_min_err_coeff.F90 diff --git a/src/All_Coeffs_mod.F90 b/src/All_Coeffs_mod.F90 index ef965d8..76dfcc9 100644 --- a/src/All_Coeffs_mod.F90 +++ b/src/All_Coeffs_mod.F90 @@ -194,6 +194,45 @@ module All_Coeffs_mod CCTK_REAL, dimension(3), intent(OUT) :: a CCTK_REAL, dimension(6,4), intent(OUT) :: q + a(1) = -2.500000000000000000000000000000000000000_wp + a(2) = 1.333333333333333333333333333333333333333_wp + a(3) = -0.08333333333333333333333333333333333333333_wp + + q(1,1) = 2.000000000000000000000000000000000000000_wp + q(2,1) = -5.000000000000000000000000000000000000000_wp + q(3,1) = 4.000000000000000000000000000000000000000_wp + q(4,1) = -1.000000000000000000000000000000000000000_wp + q(5,1) = 0_wp + q(6,1) = 0_wp + + q(1,2) = 1.000000000000000000000000000000000000000_wp + q(2,2) = -2.000000000000000000000000000000000000000_wp + q(3,2) = 1.000000000000000000000000000000000000000_wp + q(4,2) = 0_wp + q(5,2) = 0_wp + q(6,2) = 0_wp + + q(1,3) = -0.09302325581395348837209302325581395348837_wp + q(2,3) = 1.372093023255813953488372093023255813953_wp + q(3,3) = -2.558139534883720930232558139534883720930_wp + q(4,3) = 1.372093023255813953488372093023255813953_wp + q(5,3) = -0.09302325581395348837209302325581395348837_wp + q(6,3) = 0_wp + + q(1,4) = -0.02040816326530612244897959183673469387755_wp + q(2,4) = 0_wp + q(3,4) = 1.204081632653061224489795918367346938776_wp + q(4,4) = -2.408163265306122448979591836734693877551_wp + q(5,4) = 1.306122448979591836734693877551020408163_wp + q(6,4) = -0.08163265306122448979591836734693877551020_wp + + end subroutine coeffs_2_4_2 + + subroutine coeffs_2_4_2_opt ( a, q ) + + CCTK_REAL, dimension(3), intent(OUT) :: a + CCTK_REAL, dimension(6,4), intent(OUT) :: q + a(1) = -2.5_wp a(2) = 1.333333333333333333333333333333333333333_wp a(3) = -0.08333333333333333333333333333333333333333_wp @@ -226,7 +265,7 @@ module All_Coeffs_mod q(5,4) = -0.9258553965216912759079296546576971585129_wp q(6,4) = 0.4150007042443253806547133319378374518503_wp - end subroutine coeffs_2_4_2 + end subroutine coeffs_2_4_2_opt subroutine coeffs_1_6_3 ( a, q ) diff --git a/src/Derivatives2_4_2.F90 b/src/Derivatives2_4_2.F90 index 1827966..443775f 100644 --- a/src/Derivatives2_4_2.F90 +++ b/src/Derivatives2_4_2.F90 @@ -179,10 +179,10 @@ subroutine deriv2_gf_4_2 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar2 ) q(5,3) * var(:,:,nk-4) + & q(6,3) * var(:,:,nk-5) ) * idel2 dvar2(:,:,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(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) ) * idel2 dvar2(:,:,nk) = ( q(1,1) * var(:,:,nk) + & q(2,1) * var(:,:,nk-1) + & diff --git a/src/Derivatives2_4_2_min_err_coeff.F90 b/src/Derivatives2_4_2_min_err_coeff.F90 new file mode 100644 index 0000000..ed3d81f --- /dev/null +++ b/src/Derivatives2_4_2_min_err_coeff.F90 @@ -0,0 +1,177 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine deriv2_gf_4_2_opt ( var, ni, nj, nk, dir, bb, gsize, delta, dvar2 ) + + use All_Coeffs_mod + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + 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_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar2 + + CCTK_REAL, dimension(3), save :: a + CCTK_REAL, dimension(6,4), save :: q + CCTK_REAL :: idel2 + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + call coeffs_2_4_2_opt ( a, q ) + first = .false. + end if + + idel2 = 1.0_wp / delta**2 + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dvar2(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + & + q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) ) * idel2 + dvar2(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + & + q(3,2) * var(3,:,:) ) * idel2 + dvar2(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,:,:) ) * idel2 + dvar2(4,:,:) = ( q(1,4) * var(1,:,:) + q(3,4) * var(3,:,:) + & + q(4,4) * var(4,:,:) + q(5,4) * var(5,:,:) + & + q(6,4) * var(6,:,:) ) * idel2 + il = 5 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + dvar2(ni-3,:,:) = ( q(1,4) * var(ni,:,:) + & + q(3,4) * var(ni-2,:,:) + & + q(4,4) * var(ni-3,:,:) + & + q(5,4) * var(ni-4,:,:) + & + q(6,4) * var(ni-5,:,:) ) * idel2 + dvar2(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,:,:) ) * idel2 + dvar2(ni-1,:,:) = ( q(1,2) * var(ni,:,:) + & + q(2,2) * var(ni-1,:,:) + & + q(3,2) * var(ni-2,:,:) ) * idel2 + dvar2(ni,:,:) = ( q(1,1) * var(ni,:,:) + & + q(2,1) * var(ni-1,:,:) + & + q(3,1) * var(ni-2,:,:) + & + q(4,1) * var(ni-3,:,:) ) * idel2 + ir = ni - 4 + end if + dvar2(il:ir,:,:) = ( a(1) * var(il:ir,:,:) + & + a(2) * ( var(il+1:ir+1,:,:) + & + var(il-1:ir-1,:,:) ) + & + a(3) * ( var(il+2:ir+2,:,:) + & + var(il-2:ir-2,:,:) ) ) * idel2 + case (1) direction + if ( zero_derivs_y /= 0 ) then + dvar2 = zero + else + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + dvar2(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + & + q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) ) * idel2 + dvar2(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + & + q(3,2) * var(:,3,:) ) * idel2 + dvar2(:,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,:) ) * idel2 + dvar2(:,4,:) = ( q(1,4) * var(:,1,:) + q(3,4) * var(:,3,:) + & + q(4,4) * var(:,4,:) + q(5,4) * var(:,5,:) + & + q(6,4) * var(:,6,:) ) * idel2 + jl = 5 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + dvar2(:,nj-3,:) = ( q(1,4) * var(:,nj,:) + & + q(3,4) * var(:,nj-2,:) + & + q(4,4) * var(:,nj-3,:) + & + q(5,4) * var(:,nj-4,:) + & + q(6,4) * var(:,nj-5,:) ) * idel2 + dvar2(:,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,:) ) * idel2 + dvar2(:,nj-1,:) = ( q(1,2) * var(:,nj,:) + & + q(2,2) * var(:,nj-1,:) + & + q(3,2) * var(:,nj-2,:) ) * idel2 + dvar2(:,nj,:) = ( q(1,1) * var(:,nj,:) + & + q(2,1) * var(:,nj-1,:) + & + q(3,1) * var(:,nj-2,:) + & + q(4,1) * var(:,nj-3,:) ) * idel2 + jr = nj - 4 + end if + dvar2(:,jl:jr,:) = ( a(1) * var(:,jl:jr,:) + & + a(2) * ( var(:,jl+1:jr+1,:) + & + var(:,jl-1:jr-1,:) ) + & + a(3) * ( var(:,jl+2:jr+2,:) + & + var(:,jl-2:jr-2,:) ) ) * idel2 + + end if + case (2) direction + if ( zero_derivs_z /= 0 ) then + dvar2 = zero + else + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + dvar2(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + & + q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) ) * idel2 + dvar2(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + & + q(3,2) * var(:,:,3) ) * idel2 + dvar2(:,:,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) ) * idel2 + dvar2(:,:,4) = ( q(1,4) * var(:,:,1) + q(3,4) * var(:,:,3) + & + q(4,4) * var(:,:,4) + q(5,4) * var(:,:,5) + & + q(6,4) * var(:,:,6) ) * idel2 + kl = 5 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + dvar2(:,:,nk-3) = ( q(1,4) * var(:,:,nk) + & + q(3,4) * var(:,:,nk-2) + & + q(4,4) * var(:,:,nk-3) + & + q(5,4) * var(:,:,nk-4) + & + q(6,4) * var(:,:,nk-5) ) * idel2 + dvar2(:,:,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) ) * idel2 + dvar2(:,:,nk-1) = ( q(1,2) * var(:,:,nk) + & + q(2,2) * var(:,:,nk-1) + & + q(3,2) * var(:,:,nk-2) ) * idel2 + dvar2(:,:,nk) = ( q(1,1) * var(:,:,nk) + & + q(2,1) * var(:,:,nk-1) + & + q(3,1) * var(:,:,nk-2) + & + q(4,1) * var(:,:,nk-3) ) * idel2 + kr = nk - 4 + end if + dvar2(:,:,kl:kr) = ( a(1) * var(:,:,kl:kr) + & + a(2) * ( var(:,:,kl+1:kr+1) + & + var(:,:,kl-1:kr-1) ) + & + a(3) * ( var(:,:,kl+2:kr+2) + & + var(:,:,kl-2:kr-2) ) ) * idel2 + end if + end select direction +end subroutine deriv2_gf_4_2_opt diff --git a/src/call_derivs2.c b/src/call_derivs2.c index 23a5899..cf193c6 100644 --- a/src/call_derivs2.c +++ b/src/call_derivs2.c @@ -41,6 +41,15 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, const CCTK_INT *gsize, const CCTK_REAL *delta, CCTK_REAL *dvar2); + void CCTK_FCALL CCTK_FNAME(deriv2_gf_4_2_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, + CCTK_REAL *dvar2); ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2]; @@ -85,17 +94,32 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, } if ( CCTK_Equals(norm_type,"Diagonal") ) { - switch(loc_order) { - case 2: { - CCTK_FNAME(deriv2_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); - break; - } - case 4: { - CCTK_FNAME(deriv2_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); - break; - } - default: + if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) { + switch(loc_order) { + case 2: { + CCTK_FNAME(deriv2_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); + break; + } + case 4: { + CCTK_FNAME(deriv2_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); + break; + } + default: CCTK_WARN (0, "Unknown 2nd derivative stencil specified"); + } + } else { + switch(loc_order) { + case 2: { + CCTK_FNAME(deriv2_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); + break; + } + case 4: { + CCTK_FNAME(deriv2_gf_4_2_opt)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); + break; + } + default: + CCTK_WARN (0, "Unknown stencil specified"); + } } } else { switch(loc_order) { diff --git a/src/make.code.defn b/src/make.code.defn index b5bb131..4357264 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -11,6 +11,7 @@ SRCS = call_derivs.c \ Derivatives2_2_1.F90 \ Derivatives_4_2.F90 \ Derivatives2_4_2.F90 \ + Derivatives2_4_2_min_err_coeff.F90 \ Derivatives_6_3.F90 \ Derivatives_6_3_min_err_coeff.F90 \ Derivatives_8_4.F90 \ diff --git a/src/make.code.deps b/src/make.code.deps index 23655a1..6091656 100644 --- a/src/make.code.deps +++ b/src/make.code.deps @@ -5,6 +5,7 @@ Coefficients2_2_1.o: All_Coeffs_mod.F90.o Derivatives_4_2.o: All_Coeffs_mod.F90.o Coefficients_4_2.o: All_Coeffs_mod.F90.o Derivatives2_4_2.o: All_Coeffs_mod.F90.o +Derivatives2_4_2_min_err_coeff.o: All_Coeffs_mod.F90.o Coefficients2_4_2.o: All_Coeffs_mod.F90.o Derivatives_4_3.o: All_Coeffs_mod.F90.o Coefficients_4_3.o: All_Coeffs_mod.F90.o -- cgit v1.2.3