From 3d5b284d6cc60be8c2685e2eddcab6b1b4801c3b Mon Sep 17 00:00:00 2001 From: diener Date: Wed, 31 Oct 2007 17:37:57 +0000 Subject: Add 8-4 second derivative operators. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@93 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- src/All_Coeffs_mod.F90 | 148 ++++++++++++++-- src/Coefficients2_8_4.F90 | 98 +++++++++++ src/Derivatives2_8_4.F90 | 441 ++++++++++++++++++++++++++++++++++++++++++++++ src/call_derivs2.c | 17 ++ src/get_coeffs2.c | 21 ++- src/make.code.defn | 2 + src/make.code.deps | 2 + 7 files changed, 706 insertions(+), 23 deletions(-) create mode 100644 src/Coefficients2_8_4.F90 create mode 100644 src/Derivatives2_8_4.F90 diff --git a/src/All_Coeffs_mod.F90 b/src/All_Coeffs_mod.F90 index 69e58a9..1b773ad 100644 --- a/src/All_Coeffs_mod.F90 +++ b/src/All_Coeffs_mod.F90 @@ -202,25 +202,25 @@ module All_Coeffs_mod 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(5,1) = zero + q(6,1) = zero 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(4,2) = zero + q(5,2) = zero + q(6,2) = zero 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(6,3) = zero q(1,4) = -0.02040816326530612244897959183673469387755_wp - q(2,4) = 0_wp + q(2,4) = zero q(3,4) = 1.204081632653061224489795918367346938776_wp q(4,4) = -2.408163265306122448979591836734693877551_wp q(5,4) = 1.306122448979591836734693877551020408163_wp @@ -566,9 +566,9 @@ module All_Coeffs_mod 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(7,1) = zero + q(8,1) = zero + q(9,1) = zero q(1,2) = 1.053412969283276450511945392491467576792_wp q(2,2) = -2.350398179749715585893060295790671217292_wp @@ -576,9 +576,9 @@ module All_Coeffs_mod 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(7,2) = zero + q(8,2) = zero + q(9,2) = zero q(1,3) = -0.6441780400836099840157383499323742776343_wp q(2,3) = 4.137556867084716586745358416328538054838_wp @@ -586,9 +586,9 @@ module All_Coeffs_mod 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(7,3) = zero + q(8,3) = zero + q(9,3) = zero q(1,4) = 0.2133575915904708589911053057162405921503_wp q(2,4) = -1.159078186228774025004665049449524164956_wp @@ -597,8 +597,8 @@ module All_Coeffs_mod 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(8,4) = zero + q(9,4) = zero q(1,5) = -0.1790783293131903008759680081249206550717_wp q(2,5) = 0.9156510515847827006897719097795268926410_wp @@ -608,7 +608,7 @@ module All_Coeffs_mod q(6,5) = 1.735582497566755532986331513689644957894_wp q(7,5) = -0.1645296432652024882569506157166433921544_wp q(8,5) = 0.01218738098260759172273708264567728830773_wp - q(9,5) = 0_wp + q(9,5) = zero q(1,6) = 0.04002001476374207590389869333272451922711_wp q(2,6) = -0.1875223548929628699496967344733377472356_wp @@ -1140,6 +1140,116 @@ module All_Coeffs_mod end subroutine coeffs_up_8_4_opt + subroutine coeffs_2_8_4 ( a, q ) + + CCTK_REAL, dimension(5), intent(OUT) :: a + CCTK_REAL, dimension(12,8), intent(OUT) :: q + + a(1) = -2.847222222222222222222222222222222222222_wp + a(2) = 1.600000000000000000000000000000000000000_wp + a(3) = -0.2000000000000000000000000000000000000000_wp + a(4) = 0.02539682539682539682539682539682539682540_wp + a(5) = -0.001785714285714285714285714285714285714286_wp + + q(1,1) = 3.700125961665716499242827594045243158419_wp + q(2,1) = -12.53408910332763232879029889760479228385_wp + q(3,1) = 17.08522275831908082197574724401198070962_wp + q(4,1) = -12.00251923331432998485655188090486316838_wp + q(5,1) = 4.335222758319080821975747244011980709621_wp + q(6,1) = -0.5340891033276323287902988976047922838482_wp + q(7,1) = -0.04987403833428350075717240595475684158085_wp + q(8,1) = zero + q(9,1) = zero + q(10,1) = zero + q(11,1) = zero + q(12,1) = zero + q(1,2) = 0.8545561806165568782871327076593746414458_wp + q(2,2) = -1.385651335219520350522366987732872958967_wp + q(3,2) = 0.03489488503609432577108173221703695010854_wp + q(4,2) = 0.6174959481995095555971180534964638500300_wp + q(5,2) = -0.01537226034806520970159454957687817246322_wp + q(6,2) = -0.1687175231686941483830240392722911697288_wp + q(7,2) = 0.07110835640429802975122382498579196986742_wp + q(8,2) = -0.008314251520179080799570741776625110292494_wp + q(9,2) = zero + q(10,2) = zero + q(11,2) = zero + q(12,2) = zero + q(1,3) = 0.1486597947505276786174351637623847280273_wp + q(2,3) = 0.2067945264072191961776578930368936098873_wp + q(3,3) = -0.6126228482043964280320857972314409478969_wp + q(4,3) = 0.6747701953118921125519965255870017768863_wp + q(5,3) = -1.911835643616456845050516722454935307278_wp + q(6,3) = 2.589340655152612946529422151593721302602_wp + q(7,3) = -1.360526641378450595342844756571493806946_wp + q(8,3) = 0.2654199615770519345489355422778686447177_wp + q(9,3) = zero + q(10,3) = zero + q(11,3) = zero + q(12,3) = zero + q(1,4) = -0.1146188247411133220554202699779770824365_wp + q(2,4) = 0.5239525746096007180691923905405359669894_wp + q(3,4) = 0.09661320523910878808200465962763626256375_wp + q(4,4) = -1.414029112733921772841533040350356265715_wp + q(5,4) = 1.222591772291551121101947203543540982713_wp + q(6,4) = -0.6020259924428416149507501533753687531597_wp + q(7,4) = 0.3679434182813619635245551059859820833417_wp + q(8,4) = -0.08042704050374588092999589599399319429637_wp + q(9,4) = zero + q(10,4) = zero + q(11,4) = zero + q(12,4) = zero + q(1,5) = 0.06885441184905793064871518875563556079618_wp + q(2,5) = -0.05682897197659991639833322507928914093129_wp + q(3,5) = -1.192628972259900927264800832893135534458_wp + q(4,5) = 5.326668514854093024959085487041627095051_wp + q(5,5) = -8.658036133958342620625148610263805114859_wp + q(6,5) = 6.351205788124211735819368535961188726091_wp + q(7,5) = -2.217166245530019638125083881812697165497_wp + q(8,5) = 0.3822584308534509596846452244576658437999_wp + q(9,5) = -0.004326821955950548698447886167190269992355_wp + q(10,5) = zero + q(11,5) = zero + q(12,5) = zero + q(1,6) = 0.03324412543986457361212616687247594270027_wp + q(2,6) = -0.2013444667830157913459721318051963127270_wp + q(3,6) = 0.5214244967518298047750888317723818148279_wp + q(4,6) = -0.8467142789066131514432621699356515897162_wp + q(5,6) = 2.050234909410534381962583240886249492721_wp + q(6,6) = -2.978298213712682019535273900331397359519_wp + q(7,6) = 1.578055555463864252208238842945777112098_wp + q(8,6) = -0.1750701723207555800161607707603235528915_wp + q(9,6) = 0.01986478753019001522837715937418159597272_wp + q(10,6) = -0.001396742873216485445745269018497143466832_wp + q(11,6) = zero + q(12,6) = zero + q(1,7) = -0.01592923134237101475574629012747169731725_wp + q(2,7) = 0.1175046087844397827007651199218350175991_wp + q(3,7) = -0.3793709248668048299249134817817685155210_wp + q(4,7) = 0.7165680376452736733332948640656278517079_wp + q(5,7) = -0.9910611454636313458152579274571837912909_wp + q(6,7) = 2.185128799543527842172503734448943570849_wp + q(7,7) = -3.127418850946242350478043794675628822387_wp + q(8,7) = 1.685621415770389829447415862701279666019_wp + q(9,7) = -0.2166153552278720352907291696202456084323_wp + q(10,7) = 0.02750671177496787749723545011050737884854_wp + q(11,7) = -0.001934065671677428886524367585895050075288_wp + q(12,7) = zero + q(1,8) = zero + q(2,8) = -0.01256791782038527350173477151453027289777_wp + q(3,8) = 0.06770121895926724008958693099893390610910_wp + q(4,8) = -0.1432796115427160189070360250265966875352_wp + q(5,8) = 0.1563022312812120974195915804261720302791_wp + q(6,8) = -0.2217546957854982665994400260882731886207_wp + q(7,8) = 1.541934599402382617061776818198730931660_wp + q(8,8) = -2.798782027001839321767451883871457955840_wp + q(9,8) = 1.585203927110954750232022339670564355947_wp + q(10,8) = -0.1981504908888693437790027924588205444934_wp + q(11,8) = 0.02516196709699928174971464031223118025313_wp + q(12,8) = -0.001769200811507761998026810646953754861548_wp + + end subroutine coeffs_2_8_4 + subroutine coeffs_1_4_3 ( a, q ) CCTK_REAL, dimension(2), intent(OUT) :: a diff --git a/src/Coefficients2_8_4.F90 b/src/Coefficients2_8_4.F90 new file mode 100644 index 0000000..b4c26dc --- /dev/null +++ b/src/Coefficients2_8_4.F90 @@ -0,0 +1,98 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine set_coeff2_8_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(5), save :: a + CCTK_REAL, dimension(12,8), save :: q + + CCTK_INT :: i, il, ir + + logical, save :: first = .true. + + if ( first ) then + call coeffs_2_8_4 ( a, q ) + first = .false. + end if + + dd = zero + imin = 0 + imax = 0 + + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dd(1:7,1) = q(1:7,1) + imin(1) = 1; imax(1) = 7 + + dd(1:8,2) = q(1:8,2) + imin(2) = 1; imax(2) = 8 + + dd(1:8,3) = q(1:8,3) + imin(3) = 1; imax(3) = 8 + + dd(1:8,4) = q(1:8,4) + imin(4) = 1; imax(4) = 8 + + dd(1:9,5) = q(1:9,5) + imin(5) = 1; imax(5) = 9 + + dd(1:10,6) = q(1:10,6) + imin(6) = 1; imax(6) = 10 + + dd(1:11,7) = q(1:11,7) + imin(7) = 1; imax(7) = 11 + + dd(2:12,8) = q(2:12,8) + imin(8) = 2; imax(8) = 12 + + il = 9 + end if + if ( bb(2) == 0 ) then + ir = nsize - gsize + else + dd(nsize-11:nsize-1,nsize-7) = q(12:2:-1,8) + imin(nsize-7) = nsize-11; imax(nsize-7) = nsize-1 + + dd(nsize-10:nsize,nsize-6) = q(11:1:-1,7) + imin(nsize-6) = nsize-10; imax(nsize-6) = nsize + + dd(nsize-9:nsize,nsize-5) = q(10:1:-1,6) + imin(nsize-5) = nsize-9; imax(nsize-5) = nsize + + dd(nsize-8:nsize,nsize-4) = q(9:1:-1,5) + imin(nsize-4) = nsize-8; imax(nsize-4) = nsize + + dd(nsize-7:nsize,nsize-3) = q(8:1:-1,4) + imin(nsize-3) = nsize-7; imax(nsize-3) = nsize + + dd(nsize-7:nsize,nsize-2) = q(8:1:-1,3) + imin(nsize-2) = nsize-7; imax(nsize-2) = nsize + + dd(nsize-7:nsize,nsize-1) = q(8:1:-1,2) + imin(nsize-1) = nsize-7; imax(nsize-1) = nsize + + dd(nsize-6:nsize,nsize) = q(7:1:-1,1) + imin(nsize) = nsize-6; imax(nsize) = nsize + + ir = nsize - 8 + end if + do i = il, ir + dd(i-4:i-1,i) = a(5:2:-1); dd(i:i+4,i) = a(1:5) + imin(i) = i-4; imax(i) = i+4 + end do +end subroutine set_coeff2_8_4 diff --git a/src/Derivatives2_8_4.F90 b/src/Derivatives2_8_4.F90 new file mode 100644 index 0000000..e47f1ae --- /dev/null +++ b/src/Derivatives2_8_4.F90 @@ -0,0 +1,441 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine deriv2_gf_8_4 ( 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(5), save :: a + CCTK_REAL, dimension(12,8), save :: q + CCTK_REAL :: idel2 + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + call coeffs_2_8_4 ( 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(5,:,:) + q(6,1) * var(6,:,:) + & + q(7,1) * var(7,:,:) ) * 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,:,:) + & + q(7,2) * var(7,:,:) + q(8,2) * var(8,:,:) ) * 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,:,:) + & + q(7,3) * var(7,:,:) + q(8,3) * var(8,:,:) ) * 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,:,:) + q(8,4) * var(8,:,:) ) * 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,5) * var(8,:,:) + & + q(9,5) * var(9,:,:) ) * 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,:,:) + q(10,6) * var(10,:,:) ) * idel2 + dvar2(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,:,:) ) * idel2 + dvar2(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,:,:) ) * idel2 + il = 9 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + dvar2(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,:,:) ) * idel2 + dvar2(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,:,:) ) * idel2 + 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,:,:) + & + q(10,6) * var(ni-9,:,:) ) * 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,:,:) + & + q(9,5) * var(ni-8,:,:) ) * 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,:,:) + & + q(8,4) * var(ni-7,:,:) ) * 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,:,:) + & + q(7,3) * var(ni-6,:,:) + & + q(8,3) * var(ni-7,:,:) ) * 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,:,:) + & + q(7,2) * var(ni-6,:,:) + & + q(8,2) * var(ni-7,:,:)) * 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,:,:) + & + q(7,1) * var(ni-6,:,:) ) * idel2 + ir = ni - 8 + 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,:,:) ) + & + a(5) * ( var(il+4:ir+4,:,:) + & + var(il-4:ir-4,:,:) ) ) * 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(:,5,:) + q(6,1) * var(:,6,:) + & + q(7,1) * var(:,7,:) ) * 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,:) + & + q(7,2) * var(:,7,:) + q(8,2) * var(:,8,:) ) * 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,:) + & + q(7,3) * var(:,7,:) + q(8,3) * var(:,8,:) ) * 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,:) + q(8,4) * var(:,8,:) ) * 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,5) * var(:,8,:) + & + q(9,5) * var(:,9,:) ) * 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,:) + q(10,6) * var(:,10,:) ) * idel2 + dvar2(:,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,:) ) * idel2 + dvar2(:,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,:) ) * idel2 + jl = 9 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + dvar2(:,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,:) ) * idel2 + dvar2(:,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,:) ) * idel2 + 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,:) + & + q(10,6) * var(:,nj-9,:) ) * 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,:) + & + q(9,5) * var(:,nj-8,:) ) * 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,:) + & + q(8,4) * var(:,nj-7,:) ) * 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,:) + & + q(7,3) * var(:,nj-6,:) + & + q(8,3) * var(:,nj-7,:) ) * 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,:) + & + q(7,2) * var(:,nj-6,:) + & + q(8,2) * var(:,nj-7,:)) * 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,:) + & + q(7,1) * var(:,nj-6,:) ) * idel2 + jr = nj - 8 + 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,:) ) + & + a(5) * ( var(:,jl+4:jr+4,:) + & + var(:,jl-4:jr-4,:) ) ) * 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(:,:,5) + q(6,1) * var(:,:,6) + & + q(7,1) * var(:,:,7) ) * 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) + & + q(7,2) * var(:,:,7) + q(8,2) * var(:,:,8) ) * 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) + & + q(7,3) * var(:,:,7) + q(8,3) * var(:,:,8) ) * 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) + q(8,4) * var(:,:,8) ) * 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,5) * var(:,:,8) + & + q(9,5) * var(:,:,9) ) * 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) + q(10,6) * var(:,:,10) ) * idel2 + dvar2(:,:,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) ) * idel2 + dvar2(:,:,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) ) * idel2 + kl = 9 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + dvar2(:,:,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) ) * idel2 + dvar2(:,:,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) ) * idel2 + 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) + & + q(10,6) * var(:,:,nk-9) ) * 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) + & + q(9,5) * var(:,:,nk-8) ) * 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) + & + q(8,4) * var(:,:,nk-7) ) * 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) + & + q(7,3) * var(:,:,nk-6) + & + q(8,3) * var(:,:,nk-7) ) * 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) + & + q(7,2) * var(:,:,nk-6) + & + q(8,2) * var(:,:,nk-7)) * 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) + & + q(7,1) * var(:,:,nk-6) ) * idel2 + kr = nk - 8 + 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) ) + & + a(5) * ( var(:,:,kl+4:kr+4) + & + var(:,:,kl-4:kr-4) ) ) * idel2 + end if + end select direction +end subroutine deriv2_gf_8_4 diff --git a/src/call_derivs2.c b/src/call_derivs2.c index 67f33f8..a216391 100644 --- a/src/call_derivs2.c +++ b/src/call_derivs2.c @@ -59,6 +59,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_8_4)(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, @@ -133,6 +142,10 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir1, CCTK_FNAME(deriv2_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); break; } + case 8: { + CCTK_FNAME(deriv2_gf_8_4)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); + break; + } default: CCTK_WARN (0, "Unknown 2nd derivative stencil specified"); } @@ -150,6 +163,10 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir1, CCTK_FNAME(deriv2_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2); break; } + case 8: { + CCTK_FNAME(deriv2_gf_8_4)(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 9d13dec..fff110d 100644 --- a/src/get_coeffs2.c +++ b/src/get_coeffs2.c @@ -53,6 +53,13 @@ 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_8_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); ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2]; @@ -107,6 +114,10 @@ void DiffCoeff2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_FNAME(set_coeff2_6_3)(&nsize,&loc_order,bb,&gsize,imin,imax,q); break; } + case 8: { + CCTK_FNAME(set_coeff2_8_4)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } default: CCTK_WARN (0, "Unknown stencil specified"); } @@ -117,20 +128,22 @@ void DiffCoeff2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, break; } case 4: { - CCTK_FNAME(set_coeff2_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + CCTK_FNAME(set_coeff2_4_2_opt)(&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; } + case 8: { + CCTK_FNAME(set_coeff2_8_4)(&nsize,&loc_order,bb,&gsize,imin,imax,q); + break; + } default: CCTK_WARN (0, "Unknown stencil specified"); } } } else { - if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) { - CCTK_WARN (0, "Unknown stencil specified"); - } + CCTK_WARN (0, "Second derivatives not implemented for restricted full norm"); } } diff --git a/src/make.code.defn b/src/make.code.defn index 5c9b052..e2a0294 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -16,6 +16,7 @@ SRCS = call_derivs.c \ Derivatives2_6_3.F90 \ Derivatives_6_3_min_err_coeff.F90 \ Derivatives_8_4.F90 \ + Derivatives2_8_4.F90 \ Derivatives_8_4_min_err_coeff.F90 \ Derivatives_4_3.F90 \ Derivatives_4_3_min_err_coeff.F90 \ @@ -52,6 +53,7 @@ SRCS = call_derivs.c \ Coefficients2_6_3.F90 \ Coefficients_6_3_min_err_coeff.F90 \ Coefficients_8_4.F90 \ + Coefficients2_8_4.F90 \ Coefficients_8_4_min_err_coeff.F90 \ Coefficients_4_3.F90 \ Coefficients_4_3_min_err_coeff.F90 \ diff --git a/src/make.code.deps b/src/make.code.deps index 523b8bf..83781e8 100644 --- a/src/make.code.deps +++ b/src/make.code.deps @@ -26,4 +26,6 @@ Derivatives_8_4.o: All_Coeffs_mod.F90.o Coefficients_8_4.o: All_Coeffs_mod.F90.o Derivatives_8_4_min_err_coeff.o: All_Coeffs_mod.F90.o Coefficients_8_4_min_err_coeff.o: All_Coeffs_mod.F90.o +Derivatives2_8_4.o: All_Coeffs_mod.F90.o +Coefficients2_8_4.o: All_Coeffs_mod.F90.o -- cgit v1.2.3