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
|
c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyFO/src/calcrhs.F77,v 1.7 2004/05/07 22:06:08 schnetter Exp $
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
function WaveToyFO_Deriv (fm2, fm1, f0, fp1, fp2, dx, dir)
implicit none
CCTK_REAL WaveToyFO_Deriv
CCTK_REAL fm2, fm1, f0, fp1, fp2
CCTK_REAL dx
integer dir
if (dir .eq. 0) then
WaveToyFO_Deriv = (fp1 - fm1) / (2 * dx)
else if (dir .eq. -1) then
c Use first order one-sided derivatives here.
c These are less dissipative than second order one-sided derivatives.
c (Sic.)
WaveToyFO_Deriv = (f0 - fm1) / dx
c$$$ WaveToyFO_Deriv = (+ 3*f0 - 4*fm1 + fm2) / (2*dx)
else if (dir .eq. +1) then
WaveToyFO_Deriv = (fp1 - f0) / dx
c$$$ WaveToyFO_Deriv = (- 3*f0 + 4*fp1 - fp2) / (2*dx)
else
call CCTK_WARN (0, "internal error")
end if
end
#define CHOOSEDIR(i,imin,imax) ((min(i+1,imax) - (i+1)) + (max(i-1,imin) - (i-1)))
#define DIFFX(var,i,j,k) WaveToyFO_Deriv(var((i)-2,j,k), var((i)-1,j,k), var(i,j,k), var((i)+1,j,k), var((i)+2,j,k), dx(1), idir(1))
#define DIFFY(var,i,j,k) WaveToyFO_Deriv(var(i,(j)-2,k), var(i,(j)-1,k), var(i,j,k), var(i,(j)+1,k), var(i,(j)+2,k), dx(2), idir(2))
#define DIFFZ(var,i,j,k) WaveToyFO_Deriv(var(i,j,(k)-2), var(i,j,(k)-1), var(i,j,k), var(i,j,(k)+1), var(i,j,(k)+2), dx(3), idir(3))
subroutine WaveToyFO_CalcRHS (CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
CCTK_REAL dx(3)
integer bndwidth(3)
integer imin(3), imax(3), idir(3)
integer i, j, k
integer d
external WaveToyFO_Deriv
CCTK_REAL WaveToyFO_Deriv
do d=1,3
bndwidth(d) = 0
end do
do d=1,3
if (cctk_bbox(2*d).ne.0) then
imin(d) = 1+bndwidth(d)
else
imin(d) = 1+cctk_nghostzones(d)
end if
if (cctk_bbox(2*d+1).ne.0) then
imax(d) = cctk_lsh(d)-bndwidth(d)
else
imax(d) = cctk_lsh(d)-cctk_nghostzones(d)
end if
end do
do d=1,3
dx(d) = CCTK_DELTA_SPACE(3)
end do
do k=imin(3),imax(3)
idir(3) = CHOOSEDIR(k,imin(3),imax(3))
do j=imin(2),imax(2)
idir(2) = CHOOSEDIR(j,imin(2),imax(2))
do i=imin(1),imax(1)
idir(1) = CHOOSEDIR(i,imin(1),imax(1))
phidot(i,j,k) = DIFFX (psix,i,j,k)
$ + DIFFY (psiy,i,j,k)
$ + DIFFZ (psiz,i,j,k)
psixdot(i,j,k) = DIFFX (phi,i,j,k)
psiydot(i,j,k) = DIFFY (phi,i,j,k)
psizdot(i,j,k) = DIFFZ (phi,i,j,k)
end do
end do
end do
end
|
