1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
|
! $Header$
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
subroutine SBP_GetScalProdDiag ( cctkGH, dir, nsize, sigmad )
implicit none
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
CCTK_POINTER_TO_CONST, intent(IN) :: cctkGH
CCTK_INT, intent(IN) :: dir
CCTK_INT, intent(IN) :: nsize
CCTK_REAL, dimension(nsize), intent(OUT) :: sigmad
CCTK_REAL, parameter :: zero = 0.0
integer, parameter :: wp = kind(zero)
integer, dimension(6) :: onesided
CCTK_REAL, dimension(2), parameter :: bmask_2 = (/ 0.5_wp, 1.0_wp /)
CCTK_REAL, dimension(4), parameter :: bmask_4 = (/ 17.0_wp/48.0_wp, &
59.0_wp/48.0_wp, &
43.0_wp/48.0_wp, &
49.0_wp/48.0_wp /)
CCTK_REAL, dimension(6), parameter :: bmask_6 = (/ 13649.0_wp/43200._wp, &
12013.0_wp/8640._wp, &
2711.0_wp/4320.0_wp, &
5359.0_wp/4320.0_wp, &
7877.0_wp/8640.0_wp, &
43801.0_wp/43200.0_wp /)
CCTK_REAL, dimension(8), parameter :: bmask_8 = (/ 1498139.0_wp/5080320.0_wp,&
1107307.0_wp/725760.0_wp, &
20761.0_wp/80640.0_wp, &
1304999.0_wp/725760.0_wp, &
299527.0_wp/725760.0_wp, &
103097.0_wp/80640.0_wp, &
670091.0_wp/725760.0_wp, &
5127739.0_wp/5080320.0_wp/)
if ( dir < 0 .or. dir > 2 ) then
call CCTK_WARN(0, 'Error: dir is outside the legal range')
end if
call SBP_determine_onesided_stencil (cctkGH, onesided)
sigmad = 1.0_wp
if ( onesided(dir*2+1) == 1 ) then
select case (order)
case (2)
sigmad(1:2) = bmask_2
case (4)
sigmad(1:4) = bmask_4
case (6)
sigmad(1:6) = bmask_6
case (8)
sigmad(1:8) = bmask_8
end select
end if
if ( onesided(dir*2+2) == 1 ) then
select case (order)
case (2)
sigmad(nsize:nsize-1:-1) = bmask_2
case (4)
sigmad(nsize:nsize-3:-1) = bmask_4
case (6)
sigmad(nsize:nsize-5:-1) = bmask_6
case (8)
sigmad(nsize:nsize-7:-1) = bmask_8
end select
end if
end subroutine SBP_GetScalProdDiag
|
