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
|
! $Header$
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
subroutine NoExcision_Reduce (cctk_iteration, cctk_lsh, rhs, x, y, z)
implicit none
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
integer :: cctk_iteration
integer :: cctk_lsh(3)
CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: rhs, x, y, z
CCTK_REAL, parameter :: zero=0
CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: dist2
CCTK_REAL :: cx, cy, cz, radx, rady, radz, width
integer :: n
integer, parameter :: sm_linear = 1
integer, parameter :: sm_spline = 2
integer :: sm_type
do n = 1, num_regions
if (reduce_rhs(n) /= 0) then
cx = centre_x(n)
cy = centre_y(n)
cz = centre_z(n)
if (CCTK_EQUALS(region_shape(n), "sphere")) then
radx = radius(n)
rady = radius(n)
radz = radius(n)
else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
radx = radius_x(n)
rady = radius_y(n)
radz = radius_z(n)
else
call CCTK_WARN (0, "internal error")
end if
if (CCTK_EQUALS (smoothing_function(n), "linear")) then
sm_type = sm_linear
else if (CCTK_EQUALS (smoothing_function(n), "spline")) then
sm_type = sm_spline
else
call CCTK_WARN (0, "internal error")
end if
width = smoothing_zone_width(n)
dist2 = (x - cx)**2 + (y - cy)**2 + (z - cz)**2
! TODO: reduce not only the RHS, but also drive the variables
! towards Minkowski
where (dist2 <= 1)
rhs = smooth (rhs, rhs * reduction_factor(n), dist2)
end where
end if
end do
contains
elemental function smooth (oldval, goodval, dist2) result (newval)
CCTK_REAL :: newval
CCTK_REAL, intent(in) :: oldval
CCTK_REAL, intent(in) :: goodval
CCTK_REAL, intent(in) :: dist2
CCTK_REAL, parameter :: two=2, half=1/two
CCTK_REAL :: dist
CCTK_REAL :: x, f
if (dist2 <= (1 - width)**2) then
newval = goodval
else if (width == 0 .or. dist2 >= 1) then
newval = oldval
else
dist = sqrt(dist2)
x = 1 - (1 - dist) / width
select case (sm_type)
case (sm_linear)
f = x
case (sm_spline)
if (x < half) then
f = 2 * x**2
else
f = 1 - 2 * (1-x)**2
end if
end select
newval = oldval * f + goodval * (1 - f)
end if
end function smooth
end subroutine NoExcision_Reduce
|
