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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
|
! $Header$
#include "cctk.h"
#include "cctk_Functions.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
subroutine SBP_SetNormMask (CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
CCTK_REAL, parameter :: zero = 0.0
integer, parameter :: wp = kind(zero)
integer :: onesided(6)
CCTK_INT :: np
CCTK_REAL, dimension(:), allocatable :: mask_x, mask_y, mask_z
! Note: The first number is twice the value from the paper, since Carpet
! Multiplies by 1/2 at the boundary when doing the sum reduction.
CCTK_REAL, dimension(2), parameter :: bmask_2 = (/ 1.0_wp, 1.0_wp /)
CCTK_REAL, dimension(3), parameter :: bmask_3 = (/ 1.0_wp, 1.0_wp, 1.0_wp /)
CCTK_REAL, dimension(4), parameter :: bmask_4 = (/ 34.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 = (/ 27298.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 = (/ 2996278.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/)
CCTK_REAL, dimension(8) :: bmask
CCTK_INT :: ni, nj, nk
ni = cctk_lsh(1); nj = cctk_lsh(2); nk = cctk_lsh(3)
call SBP_determine_onesided_stencil (cctkGH, onesided)
if (any (onesided/=0)) then
allocate ( mask_x(ni), mask_y(nj), mask_z(nk) )
mask_x = 1.0d0; mask_y = 1.0d0; mask_z = 1.0d0
if ( CCTK_EQUALS(norm_type,'Diagonal') ) then
select case (order)
case (2)
bmask(1:2) = bmask_2
case (4)
bmask(1:4) = bmask_4
case(6)
bmask(1:6) = bmask_6
case(8)
bmask(1:8) = bmask_8
case default
call CCTK_WARN (0, "Unknown stencil specified")
end select
np = order
else
select case (order)
case (4,6)
bmask(1:3) = bmask_3
np = 3
case default
call CCTK_WARN (0, "Unknown stencil specified")
end select
end if
if (onesided(1)/=0) mask_x(1:np) = bmask(1:np)
if (onesided(2)/=0) mask_x(ni:ni-np+1:-1) = bmask(1:np)
if (onesided(3)/=0) mask_y(1:np) = bmask(1:np)
if (onesided(4)/=0) mask_y(nj:nj-np+1:-1) = bmask(1:np)
if (onesided(5)/=0) mask_z(1:np) = bmask(1:np)
if (onesided(6)/=0) mask_z(nk:nk-np+1:-1) = bmask(1:np)
nmask = outer_prod_3d ( mask_x, mask_y, mask_z )
deallocate ( mask_x, mask_y, mask_z )
else
nmask = 1.0_wp
end if
contains
function outer_prod_3d ( ax, ay, az )
CCTK_REAL, dimension(:), intent(IN) :: ax
CCTK_REAL, dimension(:), intent(IN) :: ay
CCTK_REAL, dimension(:), intent(IN) :: az
CCTK_REAL, dimension(size(ax),size(ay),size(az)) :: outer_prod_3d
CCTK_INT :: nx, ny, nz, i, j, k
nx = size(ax); ny = size(ay); nz = size(az)
do k = 1, nz
do j = 1, ny
do i = 1, nx
outer_prod_3d(i,j,k) = ax(i) * ay(j) * az(k)
end do
end do
end do
end function outer_prod_3d
end subroutine SBP_SetNormMask
|
