diff options
Diffstat (limited to 'src/Derivatives_6_3_min_err_coeff.F90')
| -rw-r--r-- | src/Derivatives_6_3_min_err_coeff.F90 | 281 |
1 files changed, 281 insertions, 0 deletions
diff --git a/src/Derivatives_6_3_min_err_coeff.F90 b/src/Derivatives_6_3_min_err_coeff.F90 new file mode 100644 index 0000000..603c66c --- /dev/null +++ b/src/Derivatives_6_3_min_err_coeff.F90 @@ -0,0 +1,281 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine deriv_gf_6_3_opt ( var, ni, nj, nk, dir, bb, gsize, delta, dvar ) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL, parameter :: zero = 0.0 + integer, parameter :: wp = kind(zero) + 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) :: dvar + + CCTK_REAL, dimension(3), save :: a + CCTK_REAL, dimension(9,6), save :: q + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + 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) = -1.582533518939116418785258993332844897062_wp + q(2,1) = 2.033426786468126253898161347360808173712_wp + q(3,1) = -0.1417052898146741610733887894481170575600_wp + q(4,1) = -0.4501096599735708523162117824920488989702_wp + q(5,1) = 0.1042956382142412661862395105494407610836_wp + q(6,1) = 0.03662604404499391209045870736276191879693_wp + q(7,1) = 0_wp + q(8,1) = 0_wp + q(9,1) = 0_wp + q(1,2) = -0.4620701275035953590186631853846278325646_wp + q(2,2) = 0_wp + q(3,2) = 0.2873679417026202568532985205129449923126_wp + q(4,2) = 0.2585974499280928196267362923074433487080_wp + q(5,2) = -0.06894808744606961472005221923058251153103_wp + q(6,2) = -0.01494717668104810274131940820517799692506_wp + q(7,2) = 0_wp + q(8,2) = 0_wp + q(9,2) = 0_wp + q(1,3) = 0.07134398748360337973038301686379010397038_wp + q(2,3) = -0.6366933020423417826592908754928085932593_wp + q(3,3) = 0_wp + q(4,3) = 0.6067199374180168986519150843189505198519_wp + q(5,3) = -0.02338660408468356531858175098561718651857_wp + q(6,3) = -0.01798401877459493040442547470431484404443_wp + q(7,3) = 0_wp + q(8,3) = 0_wp + q(9,3) = 0_wp + q(1,4) = 0.1146397975178068401430112823144985150596_wp + q(2,4) = -0.2898424301162697370942324201800071793273_wp + q(3,4) = -0.3069262456316931913128086944558079603132_wp + q(4,4) = 0_wp + q(5,4) = 0.5203848121857539166740071338174418292578_wp + q(6,4) = -0.05169127637022742348368508279860701098408_wp + q(7,4) = 0.01343534241462959507370778130248180630715_wp + q(8,4) = 0_wp + q(9,4) = 0_wp + q(1,5) = -0.03614399304268576976452921364705641609825_wp + q(2,5) = 0.1051508663818248421520867474440761344449_wp + q(3,5) = 0.01609777419666805778308369351834662756172_wp + q(4,5) = -0.7080721616106272031118456849378369336023_wp + q(5,5) = 0_wp + q(6,5) = 0.7692160858661111736140494493705980473867_wp + q(7,5) = -0.1645296432652024882569506157166433921544_wp + q(8,5) = 0.01828107147391138758410562396851593246160_wp + q(9,5) = 0_wp + q(1,6) = -0.01141318406360863692889821914555232596651_wp + q(2,6) = 0.02049729840293952857599941220163960606616_wp + q(3,6) = 0.01113095018331244864875173213474522093204_wp + q(4,6) = 0.06324365883611076515355091406993789453750_wp + q(5,6) = -0.6916640154753724474963890679085181638850_wp + q(6,6) = 0_wp + q(7,6) = 0.7397091390607520376247117645715851236273_wp + q(8,6) = -0.1479418278121504075249423529143170247255_wp + q(9,6) = 0.01643798086801671194721581699047966941394_wp + + first = .false. + end if + + idel = 1.0_wp / delta + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dvar(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,:,:) ) * idel + dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(3,2) * var(3,:,:) + & + q(4,2) * var(4,:,:) + q(5,2) * var(5,:,:) + & + q(6,2) * var(6,:,:) ) * idel + dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + & + q(4,3) * var(4,:,:) + q(5,3) * var(5,:,:) + & + q(6,3) * var(6,:,:) ) * idel + dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + & + q(3,4) * var(3,:,:) + q(5,4) * var(5,:,:) + & + q(6,4) * var(6,:,:) + q(7,4) * var(7,:,:) ) * idel + dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + & + q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + & + q(6,5) * var(6,:,:) + q(7,5) * var(7,:,:) + & + q(8,5) * var(8,:,:) ) * idel + dvar(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(7,6) * var(7,:,:) + & + q(8,6) * var(8,:,:) + q(9,6) * var(9,:,:) ) * idel + + il = 7 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + dvar(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(7,6) * var(ni-6,:,:) + & + q(8,6) * var(ni-7,:,:) + & + q(9,6) * var(ni-8,:,:) ) * idel + dvar(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(6,5) * var(ni-5,:,:) + q(7,5) * var(ni-6,:,:) + & + q(8,5) * var(ni-7,:,:) ) * idel + dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + & + q(3,4) * var(ni-2,:,:) + q(5,4) * var(ni-4,:,:) + & + q(6,4) * var(ni-5,:,:) + & + q(7,4) * var(ni-6,:,:) ) * idel + dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + & + q(4,3) * var(ni-3,:,:) + q(5,3) * var(ni-4,:,:) + & + q(6,3) * var(ni-5,:,:) ) * idel + dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(3,2) * var(ni-2,:,:) + & + q(4,2) * var(ni-3,:,:) + q(5,2) * var(ni-4,:,:) + & + q(6,2) * var(ni-5,:,:) ) * idel + dvar(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,:,:) ) * idel + + ir = ni - 6 + end if + dvar(il:ir,:,:) = ( a(1) * ( var(il+1:ir+1,:,:) - & + var(il-1:ir-1,:,:) ) + & + a(2) * ( var(il+2:ir+2,:,:) - & + var(il-2:ir-2,:,:) ) + & + a(3) * ( var(il+3:ir+3,:,:) - & + var(il-3:ir-3,:,:) ) ) * idel + case (1) direction + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + dvar(:,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,:) ) * idel + dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(3,2) * var(:,3,:) + & + q(4,2) * var(:,4,:) + q(5,2) * var(:,5,:) + & + q(6,2) * var(:,6,:) ) * idel + dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + & + q(4,3) * var(:,4,:) + q(5,3) * var(:,5,:) + & + q(6,3) * var(:,6,:) ) * idel + dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + & + q(3,4) * var(:,3,:) + q(5,4) * var(:,5,:) + & + q(6,4) * var(:,6,:) + q(7,4) * var(:,7,:) ) * idel + dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + & + q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + & + q(6,5) * var(:,6,:) + q(7,5) * var(:,7,:) + & + q(8,5) * var(:,8,:) ) * idel + dvar(:,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(7,6) * var(:,7,:) + & + q(8,6) * var(:,8,:) + q(9,6) * var(:,9,:) ) * idel + + jl = 7 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + dvar(:,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(7,6) * var(:,nj-6,:) + & + q(8,6) * var(:,nj-7,:) + & + q(9,6) * var(:,nj-8,:) ) * idel + dvar(:,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(6,5) * var(:,nj-5,:) + q(7,5) * var(:,nj-6,:) + & + q(8,5) * var(:,nj-7,:) ) * idel + dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + & + q(3,4) * var(:,nj-2,:) + q(5,4) * var(:,nj-4,:) + & + q(6,4) * var(:,nj-5,:) + & + q(7,4) * var(:,nj-6,:) ) * idel + dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + & + q(4,3) * var(:,nj-3,:) + q(5,3) * var(:,nj-4,:) + & + q(6,3) * var(:,nj-5,:) ) * idel + dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(3,2) * var(:,nj-2,:) + & + q(4,2) * var(:,nj-3,:) + q(5,2) * var(:,nj-4,:) + & + q(6,2) * var(:,nj-5,:) ) * idel + dvar(:,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,:) ) * idel + + jr = nj - 6 + end if + dvar(:,jl:jr,:) = ( a(1) * ( var(:,jl+1:jr+1,:) - & + var(:,jl-1:jr-1,:) ) + & + a(2) * ( var(:,jl+2:jr+2,:) - & + var(:,jl-2:jr-2,:) ) + & + a(3) * ( var(:,jl+3:jr+3,:) - & + var(:,jl-3:jr-3,:) ) ) * idel + case (2) direction + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + dvar(:,:,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) ) * idel + dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(3,2) * var(:,:,3) + & + q(4,2) * var(:,:,4) + q(5,2) * var(:,:,5) + & + q(6,2) * var(:,:,6) ) * idel + dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + & + q(4,3) * var(:,:,4) + q(5,3) * var(:,:,5) + & + q(6,3) * var(:,:,6) ) * idel + dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + & + q(3,4) * var(:,:,3) + q(5,4) * var(:,:,5) + & + q(6,4) * var(:,:,6) + q(7,4) * var(:,:,7) ) * idel + dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + & + q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + & + q(6,5) * var(:,:,6) + q(7,5) * var(:,:,7) + & + q(8,5) * var(:,:,8) ) * idel + dvar(:,:,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(7,6) * var(:,:,7) + & + q(8,6) * var(:,:,8) + q(9,6) * var(:,:,9) ) * idel + + kl = 7 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + dvar(:,:,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(7,6) * var(:,:,nk-6) + & + q(8,6) * var(:,:,nk-7) + & + q(9,6) * var(:,:,nk-8) ) * idel + dvar(:,:,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(6,5) * var(:,:,nk-5) + q(7,5) * var(:,:,nk-6) + & + q(8,5) * var(:,:,nk-7) ) * idel + dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + & + q(3,4) * var(:,:,nk-2) + q(5,4) * var(:,:,nk-4) + & + q(6,4) * var(:,:,nk-5) + & + q(7,4) * var(:,:,nk-6) ) * idel + dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + & + q(4,3) * var(:,:,nk-3) + q(5,3) * var(:,:,nk-4) + & + q(6,3) * var(:,:,nk-5) ) * idel + dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(3,2) * var(:,:,nk-2) + & + q(4,2) * var(:,:,nk-3) + q(5,2) * var(:,:,nk-4) + & + q(6,2) * var(:,:,nk-5) ) * idel + dvar(:,:,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) ) * idel + + kr = nk - 6 + end if + dvar(:,:,kl:kr) = ( a(1) * ( var(:,:,kl+1:kr+1) - & + var(:,:,kl-1:kr-1) ) + & + a(2) * ( var(:,:,kl+2:kr+2) - & + var(:,:,kl-2:kr-2) ) + & + a(3) * ( var(:,:,kl+3:kr+3) - & + var(:,:,kl-3:kr-3) ) ) * idel + end select direction +end subroutine deriv_gf_6_3_opt |