From b05ae4f1041d3c1c3381a093170c1fd805fc4517 Mon Sep 17 00:00:00 2001 From: diener Date: Fri, 19 Oct 2007 09:24:47 +0000 Subject: Added 6-3 minimal bandwidth 2nd derivative operators. Added 4-2 minimal bandwidth 2nd derivative operators. Moved the 4-2 optimised operators. Fixed missing break statements get_up_coeffs.c. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@89 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- src/All_Coeffs_mod.F90 | 77 +++++++++ src/Coefficients2_4_2.F90 | 24 +-- src/Coefficients2_4_2_min_err_coeff.F90 | 74 ++++++++ src/Coefficients2_6_3.F90 | 86 ++++++++++ src/Derivatives2_6_3.F90 | 294 ++++++++++++++++++++++++++++++++ src/call_derivs2.c | 17 ++ src/get_coeffs2.c | 22 +++ src/get_up_coeffs.c | 2 + src/make.code.defn | 3 + src/make.code.deps | 5 +- 10 files changed, 591 insertions(+), 13 deletions(-) create mode 100644 src/Coefficients2_4_2_min_err_coeff.F90 create mode 100644 src/Coefficients2_6_3.F90 create mode 100644 src/Derivatives2_6_3.F90 diff --git a/src/All_Coeffs_mod.F90 b/src/All_Coeffs_mod.F90 index 76dfcc9..69e58a9 100644 --- a/src/All_Coeffs_mod.F90 +++ b/src/All_Coeffs_mod.F90 @@ -497,6 +497,7 @@ module All_Coeffs_mod q2(7,1) = zero q2(8,1) = zero q2(9,1) = zero + q2(1,2) = -0.4261090852993166609415800254745304380753_wp q2(2,2) = -0.1198701406809289935902771997003246482977_wp q2(3,2) = 0.4312121105197350491616311601533345702699_wp @@ -506,6 +507,7 @@ module All_Coeffs_mod q2(7,2) = zero q2(8,2) = zero q2(9,2) = zero + q2(1,3) = -0.008331409049041400793408941823041323547141_wp q2(2,3) = -0.3179917159117626605641230407454828831892_wp q2(3,3) = -0.5046108447067502766506824050165990409443_wp @@ -515,6 +517,7 @@ module All_Coeffs_mod q2(7,3) = zero q2(8,3) = zero q2(9,3) = zero + q2(1,4) = 0.1280751399324364352167190636169803213668_wp q2(2,4) = -0.3704544846040473075364791079948980171701_wp q2(3,4) = -0.1053961094122492652071919749185808657060_wp @@ -524,6 +527,7 @@ module All_Coeffs_mod q2(7,4) = 0.02687068482925919014741556260496361261429_wp q2(8,4) = zero q2(9,4) = zero + q2(1,5) = -0.03614399304268576976452921364705641609825_wp q2(2,5) = 0.1234319378557362297361923714125920669065_wp q2(3,5) = -0.09358865464680026772155005029274896720786_wp @@ -533,6 +537,7 @@ module All_Coeffs_mod q2(7,5) = -0.2742160721086708137615843595277389869240_wp q2(8,5) = 0.03656214294782277516821124793703186492319_wp q2(9,5) = zero + q2(1,6) = -0.01141318406360863692889821914555232596651_wp q2(2,6) = 0.02049729840293952857599941220163960606616_wp q2(3,6) = 0.02756893105132916059596754912522489034598_wp @@ -545,6 +550,78 @@ module All_Coeffs_mod end subroutine coeffs_up_6_3_opt + subroutine coeffs_2_6_3 ( a, q ) + + CCTK_REAL, dimension(4), intent(OUT) :: a + CCTK_REAL, dimension(9,6), intent(OUT) :: q + + a(1) = -2.722222222222222222222222222222222222222_wp + a(2) = 1.500000000000000000000000000000000000000_wp + a(3) = -0.1500000000000000000000000000000000000000_wp + a(4) = 0.01111111111111111111111111111111111111111_wp + + q(1,1) = 2.788238454587637677973966346740908979901_wp + q(2,1) = -8.024525606271521723203165067037878232838_wp + q(3,1) = 8.215717879209710113072996800742423132342_wp + q(4,1) = -3.382384545876376779739663467409089799008_wp + q(5,1) = 0.2745256062715217232031650670378782328376_wp + q(6,1) = 0.1284282120790289886927003199257576867658_wp + q(7,1) = 0_wp + q(8,1) = 0_wp + q(9,1) = 0_wp + + q(1,2) = 1.053412969283276450511945392491467576792_wp + q(2,2) = -2.350398179749715585893060295790671217292_wp + q(3,2) = 1.867463026166097838452787258248009101251_wp + q(4,2) = -1.034129692832764505119453924914675767918_wp + q(5,2) = 0.6003981797497155858930602957906712172924_wp + q(6,2) = -0.1367463026166097838452787258248009101251_wp + q(7,2) = 0_wp + q(8,2) = 0_wp + q(9,2) = 0_wp + + q(1,3) = -0.6441780400836099840157383499323742776343_wp + q(2,3) = 4.137556867084716586745358416328538054838_wp + q(3,3) = -8.108447067502766506824050165990409443010_wp + q(4,3) = 6.941780400836099840157383499323742776343_wp + q(5,3) = -2.887556867084716586745358416328538054838_wp + q(6,3) = 0.5608447067502766506824050165990409443010_wp + q(7,3) = 0_wp + q(8,3) = 0_wp + q(9,3) = 0_wp + + q(1,4) = 0.2133575915904708589911053057162405921503_wp + q(2,4) = -1.159078186228774025004665049449524164956_wp + q(3,4) = 3.511693723953473906823412328170678609193_wp + q(4,4) = -4.723144865335572557069104932512284630217_wp + q(5,4) = 2.489690240716551595446911737264415002799_wp + q(6,4) = -0.3414753996392361759034645767245132798408_wp + q(7,4) = 0.008956894943086396715805187534987870871431_wp + q(8,4) = 0_wp + q(9,4) = 0_wp + + q(1,5) = -0.1790783293131903008759680081249206550717_wp + q(2,5) = 0.9156510515847827006897719097795268926410_wp + q(3,5) = -1.987601032541999915365409843002835258770_wp + q(4,5) = 3.387647581566586263806017519360162498413_wp + q(5,5) = -3.719859506580339384706529558630612331260_wp + q(6,5) = 1.735582497566755532986331513689644957894_wp + q(7,5) = -0.1645296432652024882569506157166433921544_wp + q(8,5) = 0.01218738098260759172273708264567728830773_wp + q(9,5) = 0_wp + + q(1,6) = 0.04002001476374207590389869333272451922711_wp + q(2,6) = -0.1875223548929628699496967344733377472356_wp + q(3,6) = 0.3471267779274445788908929019885390744504_wp + q(4,6) = -0.4177910702190969764769449708149737829426_wp + q(5,6) = 1.560601736642238000654475164189554271972_wp + q(6,6) = -2.684870208442729618045250108445012670944_wp + q(7,6) = 1.479418278121504075249423529143170247255_wp + q(8,6) = -0.1479418278121504075249423529143170247255_wp + q(9,6) = 0.01095865391201114129814387799365311294263_wp + + end subroutine coeffs_2_6_3 + subroutine coeffs_1_8_4 ( a, q ) CCTK_REAL, dimension(4), intent(OUT) :: a diff --git a/src/Coefficients2_4_2.F90 b/src/Coefficients2_4_2.F90 index 14a1cec..5883b02 100644 --- a/src/Coefficients2_4_2.F90 +++ b/src/Coefficients2_4_2.F90 @@ -36,14 +36,14 @@ subroutine set_coeff2_4_2 ( nsize, loc_order, bb, gsize, imin, imax, dd ) if ( bb(1) == 0 ) then il = 1 + gsize else - dd(1:6,1) = q(1:6,1) - imin(1) = 1; imax(1) = 6 + dd(1:4,1) = q(1:4,1) + imin(1) = 1; imax(1) = 4 - dd(1:6,2) = q(1:6,2) - imin(2) = 1; imax(2) = 6 + dd(1:3,2) = q(1:3,2) + imin(2) = 1; imax(2) = 3 - dd(1:6,3) = q(1:6,3) - imin(3) = 1; imax(3) = 6 + dd(1:5,3) = q(1:5,3) + imin(3) = 1; imax(3) = 5 dd(1:6,4) = q(1:6,4) imin(4) = 1; imax(4) = 6 @@ -56,14 +56,14 @@ subroutine set_coeff2_4_2 ( nsize, loc_order, bb, gsize, imin, imax, dd ) dd(nsize-5:nsize,nsize-3) = q(6:1:-1,4) imin(nsize-3) = nsize-5; imax(nsize-3) = nsize - dd(nsize-5:nsize,nsize-2) = q(6:1:-1,3) - imin(nsize-2) = nsize-5; imax(nsize-2) = nsize + dd(nsize-4:nsize,nsize-2) = q(5:1:-1,3) + imin(nsize-2) = nsize-4; imax(nsize-2) = nsize - dd(nsize-5:nsize,nsize-1) = q(6:1:-1,2) - imin(nsize-1) = nsize-5; imax(nsize-1) = nsize + dd(nsize-2:nsize,nsize-1) = q(3:1:-1,2) + imin(nsize-1) = nsize-2; imax(nsize-1) = nsize - dd(nsize-5:nsize,nsize) = q(6:1:-1,1) - imin(nsize) = nsize-5; imax(nsize) = nsize + dd(nsize-3:nsize,nsize) = q(4:1:-1,1) + imin(nsize) = nsize-3; imax(nsize) = nsize ir = nsize - 4 end if diff --git a/src/Coefficients2_4_2_min_err_coeff.F90 b/src/Coefficients2_4_2_min_err_coeff.F90 new file mode 100644 index 0000000..9720d05 --- /dev/null +++ b/src/Coefficients2_4_2_min_err_coeff.F90 @@ -0,0 +1,74 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine set_coeff2_4_2_opt ( 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(6,4), save :: q + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( first ) then + call coeffs_2_4_2_opt ( a, q ) + first = .false. + end if + + dd = zero + imin = 0 + imax = 0 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dd(1:6,1) = q(1:6,1) + imin(1) = 1; imax(1) = 6 + + dd(1:6,2) = q(1:6,2) + imin(2) = 1; imax(2) = 6 + + dd(1:6,3) = q(1:6,3) + imin(3) = 1; imax(3) = 6 + + dd(1:6,4) = q(1:6,4) + imin(4) = 1; imax(4) = 6 + + il = 5 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + dd(nsize-5:nsize,nsize-3) = q(6:1:-1,4) + imin(nsize-3) = nsize-5; imax(nsize-3) = nsize + + dd(nsize-5:nsize,nsize-2) = q(6:1:-1,3) + imin(nsize-2) = nsize-5; imax(nsize-2) = nsize + + dd(nsize-5:nsize,nsize-1) = q(6:1:-1,2) + imin(nsize-1) = nsize-5; imax(nsize-1) = nsize + + dd(nsize-5:nsize,nsize) = q(6:1:-1,1) + imin(nsize) = nsize-5; imax(nsize) = nsize + + ir = nsize - 4 + end if + do i = il, ir + dd(i-2:i-1,i) = a(3:2:-1); dd(i:i+2,i) = a(1:3) + imin(i) = i-2; imax(i) = i+2 + end do +end subroutine set_coeff2_4_2_opt diff --git a/src/Coefficients2_6_3.F90 b/src/Coefficients2_6_3.F90 new file mode 100644 index 0000000..3563c48 --- /dev/null +++ b/src/Coefficients2_6_3.F90 @@ -0,0 +1,86 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine set_coeff2_6_3 ( 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(9,6), save :: q + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( first ) then + call coeffs_2_6_3 ( a, q ) + first = .false. + end if + + dd = zero + imin = 0 + imax = 0 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dd(1:6,1) = q(1:6,1) + imin(1) = 1; imax(1) = 6 + + dd(1:6,2) = q(1:6,2) + imin(2) = 1; imax(2) = 6 + + dd(1:6,3) = q(1:6,3) + imin(3) = 1; imax(3) = 6 + + dd(1:7,4) = q(1:7,4) + imin(4) = 1; imax(4) = 7 + + dd(1:8,5) = q(1:8,5) + imin(5) = 1; imax(5) = 8 + + dd(1:9,6) = q(1:9,6) + imin(6) = 1; imax(6) = 9 + + il = 7 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + dd(nsize-8:nsize,nsize-5) = q(9:1:-1,6) + imin(nsize-5) = nsize-8; imax(nsize-5) = nsize + + dd(nsize-7:nsize,nsize-4) = q(8:1:-1,5) + imin(nsize-4) = nsize-7; imax(nsize-4) = nsize + + dd(nsize-6:nsize,nsize-3) = q(7:1:-1,4) + imin(nsize-3) = nsize-6; imax(nsize-3) = nsize + + dd(nsize-5:nsize,nsize-2) = q(6:1:-1,3) + imin(nsize-2) = nsize-5; imax(nsize-2) = nsize + + dd(nsize-5:nsize,nsize-1) = q(6:1:-1,2) + imin(nsize-1) = nsize-5; imax(nsize-1) = nsize + + dd(nsize-5:nsize,nsize) = q(6:1:-1,1) + imin(nsize) = nsize-5; imax(nsize) = nsize + + ir = nsize - 6 + end if + do i = il, ir + dd(i-3:i-1,i) = a(4:2:-1); dd(i:i+3,i) = a(1:4) + imin(i) = i-3; imax(i) = i+3 + end do +end subroutine set_coeff2_6_3 diff --git a/src/Derivatives2_6_3.F90 b/src/Derivatives2_6_3.F90 new file mode 100644 index 0000000..63bf6b8 --- /dev/null +++ b/src/Derivatives2_6_3.F90 @@ -0,0 +1,294 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine deriv2_gf_6_3 ( 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(4), save :: a + CCTK_REAL, dimension(9,6), save :: q + CCTK_REAL :: idel2 + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + call coeffs_2_6_3 ( 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,:,:) + & + q(5,1) * var(4,:,:) + q(6,1) * var(6,:,:) ) * idel2 + dvar2(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,:,:) ) * 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,:,:) + q(6,3) * var(6,:,:) ) * idel2 + dvar2(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,:,:) ) * idel2 + dvar2(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,6) * var(8,:,:) ) * idel2 + dvar2(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,:,:) ) * idel2 + il = 7 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + dvar2(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,:,:) ) * idel2 + dvar2(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,:,:) ) * idel2 + dvar2(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,:,:) ) * 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,:,:) + & + q(6,3) * var(ni-5,:,:) ) * idel2 + dvar2(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,:,:) ) * 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,:,:) + & + q(5,1) * var(ni-4,:,:) + & + q(6,1) * var(ni-5,:,:) ) * idel2 + ir = ni - 6 + 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,:,:) ) + & + a(4) * ( var(il+3:ir+3,:,:) + & + var(il-3:ir-3,:,:) ) ) * 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,:) + & + q(5,1) * var(:,4,:) + q(6,1) * var(:,6,:) ) * idel2 + dvar2(:,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,:) ) * 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,:) + q(6,3) * var(:,6,:) ) * idel2 + dvar2(:,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,:) ) * idel2 + dvar2(:,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,6) * var(:,8,:) ) * idel2 + dvar2(:,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,:) ) * idel2 + jl = 7 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + dvar2(:,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,:) ) * idel2 + dvar2(:,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,:) ) * idel2 + dvar2(:,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,:) ) * 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,:) + & + q(6,3) * var(:,nj-5,:) ) * idel2 + dvar2(:,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,:) ) * 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,:) + & + q(5,1) * var(:,nj-4,:) + & + q(6,1) * var(:,nj-5,:) ) * idel2 + jr = nj - 6 + 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,:) ) + & + a(4) * ( var(:,jl+3:jr+3,:) + & + var(:,jl-3:jr-3,:) ) ) * 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) + & + q(5,1) * var(:,:,4) + q(6,1) * var(:,:,6) ) * idel2 + dvar2(:,:,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) ) * 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) + q(6,3) * var(:,:,6) ) * idel2 + dvar2(:,:,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) ) * idel2 + dvar2(:,:,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,6) * var(:,:,8) ) * idel2 + dvar2(:,:,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) ) * idel2 + kl = 7 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + dvar2(:,:,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) ) * idel2 + dvar2(:,:,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) ) * idel2 + dvar2(:,:,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) ) * 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) + & + 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(6,2) * var(:,:,nk-5) ) * 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) + & + q(5,1) * var(:,:,nk-4) + & + q(6,1) * var(:,:,nk-5) ) * idel2 + kr = nk - 6 + 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) ) + & + a(4) * ( var(:,:,kl+3:kr+3) + & + var(:,:,kl-3:kr-3) ) ) * idel2 + end if + end select direction +end subroutine deriv2_gf_6_3 diff --git a/src/call_derivs2.c b/src/call_derivs2.c index 5ec2100..67f33f8 100644 --- a/src/call_derivs2.c +++ b/src/call_derivs2.c @@ -50,6 +50,15 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir1, const CCTK_INT *gsize, const CCTK_REAL *delta, CCTK_REAL *dvar2); + void CCTK_FCALL CCTK_FNAME(deriv2_gf_6_3)(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); void CCTK_FCALL CCTK_FNAME(deriv2_mixed)(const CCTK_POINTER_TO_CONST *cctkGH, const CCTK_INT *dir1, const CCTK_INT *dir2, @@ -120,6 +129,10 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir1, CCTK_FNAME(deriv2_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); break; } + case 6: { + CCTK_FNAME(deriv2_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); + break; + } default: CCTK_WARN (0, "Unknown 2nd derivative stencil specified"); } @@ -133,6 +146,10 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir1, CCTK_FNAME(deriv2_gf_4_2_opt)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); break; } + case 6: { + CCTK_FNAME(deriv2_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); + break; + } default: CCTK_WARN (0, "Unknown stencil specified"); } diff --git a/src/get_coeffs2.c b/src/get_coeffs2.c index 9d2a0ea..9d13dec 100644 --- a/src/get_coeffs2.c +++ b/src/get_coeffs2.c @@ -39,6 +39,20 @@ void DiffCoeff2 ( 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_coeff2_4_2_opt)(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_coeff2_6_3)(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]; @@ -89,6 +103,10 @@ void DiffCoeff2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_FNAME(set_coeff2_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q); break; } + case 6: { + CCTK_FNAME(set_coeff2_6_3)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } default: CCTK_WARN (0, "Unknown stencil specified"); } @@ -102,6 +120,10 @@ void DiffCoeff2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_FNAME(set_coeff2_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q); break; } + case 6: { + CCTK_FNAME(set_coeff2_6_3)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } default: CCTK_WARN (0, "Unknown stencil specified"); } diff --git a/src/get_up_coeffs.c b/src/get_up_coeffs.c index 77226cc..3a78f6a 100644 --- a/src/get_up_coeffs.c +++ b/src/get_up_coeffs.c @@ -124,9 +124,11 @@ void DiffUpCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, } 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"); diff --git a/src/make.code.defn b/src/make.code.defn index 7c20a83..5c9b052 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -13,6 +13,7 @@ SRCS = call_derivs.c \ Derivatives2_4_2.F90 \ Derivatives2_4_2_min_err_coeff.F90 \ Derivatives_6_3.F90 \ + Derivatives2_6_3.F90 \ Derivatives_6_3_min_err_coeff.F90 \ Derivatives_8_4.F90 \ Derivatives_8_4_min_err_coeff.F90 \ @@ -46,7 +47,9 @@ SRCS = call_derivs.c \ Coefficients2_2_1.F90 \ Coefficients_4_2.F90 \ Coefficients2_4_2.F90 \ + Coefficients2_4_2_min_err_coeff.F90 \ Coefficients_6_3.F90 \ + Coefficients2_6_3.F90 \ Coefficients_6_3_min_err_coeff.F90 \ Coefficients_8_4.F90 \ Coefficients_8_4_min_err_coeff.F90 \ diff --git a/src/make.code.deps b/src/make.code.deps index 6091656..523b8bf 100644 --- a/src/make.code.deps +++ b/src/make.code.deps @@ -5,8 +5,9 @@ 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 +Derivatives2_4_2_min_err_coeff.o: All_Coeffs_mod.F90.o +Coefficients2_4_2_min_err_coeff.o: All_Coeffs_mod.F90.o Derivatives_4_3.o: All_Coeffs_mod.F90.o Coefficients_4_3.o: All_Coeffs_mod.F90.o Derivatives_4_3_min_err_coeff.o: All_Coeffs_mod.F90.o @@ -15,6 +16,8 @@ Derivatives_6_3.o: All_Coeffs_mod.F90.o Coefficients_6_3.o: All_Coeffs_mod.F90.o Derivatives_6_3_min_err_coeff.o: All_Coeffs_mod.F90.o Coefficients_6_3_min_err_coeff.o: All_Coeffs_mod.F90.o +Derivatives2_6_3.o: All_Coeffs_mod.F90.o +Coefficients2_6_3.o: All_Coeffs_mod.F90.o Derivatives_6_5.o: All_Coeffs_mod.F90.o Coefficients_6_5.o: All_Coeffs_mod.F90.o Derivatives_6_5_min_err_coeff.o: All_Coeffs_mod.F90.o -- cgit v1.2.3