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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
/*@@
@file GRHydro_MP5Reconstruct.F90
@date Fri Jan 3 2013
@author Ian Hawke, Christian Reisswig
@desc
Routines to set up the coefficient array and to perform one dimensional
M5 reconstruction of arbitrary order.
@enddesc
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
/*@@
@routine GRHydro_MP5Reconstruct1d
@date Fri Jan 3 2013
@author Christian Reisswig
@desc
Perform a one dimensional reconstruction of a given array using M5.
@enddesc
@calls
@calledby
@history
@endhistory
@@*/
#define SpaceMask_CheckStateBitsF90_1D(mask,i,type_bits,state_bits) \
(iand(mask((i)),(type_bits)).eq.(state_bits))
subroutine GRHydro_MP5Reconstruct1d(nx, v, vminus, vplus, trivial_rp, &
hydro_excision_mask)
implicit none
DECLARE_CCTK_PARAMETERS
CCTK_INT :: nx, i
CCTK_REAL, dimension(nx) :: v, vplus, vminus
CCTK_INT, dimension(nx) :: hydro_excision_mask
logical, dimension(nx) :: trivial_rp
logical, dimension(nx) :: excise
logical :: normal_m5
CCTK_REAL :: vl, vmp, djm1, dj, djp1, dm4jph, dm4jmh, vul, vav, vmd, vlc, vmin, vmax, vnorm
! sign requires its arguments to be of identical KIND
CCTK_REAL, parameter :: one = 1d0
vminus = 0.d0
vplus = 0.d0
excise = .false.
trivial_rp = .false.
!!$ Initialize excision
do i = 1, nx
if (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i) .ne. 0)) then
trivial_rp(i) = .true.
excise(i) = .true.
if (i > 1) then
trivial_rp(i-1) = .true.
end if
end if
end do
do i = 3, nx-2
!!$ Handle excision
normal_m5 = .true.
if (i < nx) then
if (excise(i+1)) then
vminus(i) = v(i)
vplus(i) = v(i)
normal_m5 = .false.
end if
end if
if (i > 1) then
if (excise(i-1)) then
vminus(i) = v(i)
vplus(i) = v(i)
normal_m5 = .false.
end if
end if
if (normal_m5) then
vnorm = sqrt(v(i-2)**2 + v(i-1)**2 + v(i)**2 + v(i+1)**2 + v(i+2)**2)
#define MINMOD(x,y) \
0.5d0*(sign(one,x) + sign(one,y)) * min(abs(x), abs(y))
#define MINMOD4(w,x,y,z) \
0.125d0*( sign(one,w)+sign(one,x) )*abs( (sign(one,w)+sign(one,y)) * (sign(one,w)+sign(one,z)) )*min(abs(w), abs(x), abs(y), abs(z))
#define MP5(am2, am1, a, ap1, ap2, arecon) \
vl = (2.0d0*am2 - 13.0d0*am1 + 47.0d0*a + 27.0d0*ap1 - 3.0d0*ap2)/60.0d0 &&\
vmp = a + MINMOD( ap1-a, mp5_alpha*(a-am1) ) &&\
if ((vl-a)*(vl-vmp) .le. mp5_eps*vnorm) then &&\
arecon = vl &&\
else &&\
djm1 = am2 -2.0d0*am1 + a &&\
dj = am1 -2.0d0*a + ap1 &&\
djp1 = a -2.d0*ap1 + ap2 &&\
dm4jph = MINMOD4(4.0d0*dj-djp1, 4.0d0*djp1-dj, dj, djp1) &&\
dm4jmh = MINMOD4(4.0d0*dj-djm1, 4.0d0*djm1-dj, dj, djm1) &&\
vul = a + mp5_alpha*(a-am1) &&\
vav = 0.5d0*(a+ap1) &&\
vmd = vav - 0.5d0*dm4jph &&\
vlc = a + 0.5d0*(a-am1) + 4.0d0/3.0d0*dm4jmh &&\
vmin = max(min(a,ap1,vmd), min(a,vul,vlc)) &&\
vmax = min(max(a,ap1,vmd), max(a,vul,vlc)) &&\
arecon = vl + MINMOD(vmin-vl, vmax-vl) &&\
endif
MP5(v(i-2), v(i-1), v(i), v(i+1), v(i+2), vplus(i))
MP5(v(i+2), v(i+1), v(i), v(i-1), v(i-2), vminus(i))
end if
end do
do i = 1, nx
if (excise(i)) then
if (i > 1) then
if (.not. excise(i-1)) then
vminus(i) = vplus(i-1)
end if
end if
vplus(i) = vminus(i)
end if
end do
do i = nx, 1, -1
if (excise(i)) then
if (i < nx) then
if (.not. excise(i+1)) then
vplus(i) = vminus(i+1)
end if
end if
vminus(i) = vplus(i)
end if
end do
end subroutine GRHydro_MP5Reconstruct1d
|
