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
|
#include <math.h>
#include "cctk.h"
#include "cctk_Parameters.h"
void BndApplyRobin3Di(cGH *GH,
int gdim,
int *sw,
int *doBC,
int *lsh,
CCTK_REAL *var,
CCTK_REAL *r,
CCTK_REAL finf,
CCTK_INT npow)
{
int xgp1,xgp2,xgp3,xgp4;
int ygp1,ygp2,ygp3,ygp4;
int zgp1,zgp2,zgp3,zgp4;
int i,j,k;
if ((sw[0]!=1) || (sw[1]!=1) || (sw[2]!=1))
CCTK_WARN(0,"Robin BC only with stencil_size=1");
/* Lower x-bound: x(2,:,:) --> x[xgp2] */
if (doBC[0]==1) {
for (k=0;k<lsh[2];k++) {
for (j=0;j<lsh[1];j++) {
xgp1 = CCTK_GFINDEX3D(GH,0,j,k);
xgp2 = CCTK_GFINDEX3D(GH,1,j,k);
xgp3 = CCTK_GFINDEX3D(GH,2,j,k);
xgp4 = CCTK_GFINDEX3D(GH,3,j,k);
var[xgp1] = finf + (
pow(r[xgp4],npow)*( var[xgp4] - finf )
-3.0*pow(r[xgp3],npow)*( var[xgp3] - finf )
+3.0*pow(r[xgp2],npow)*( var[xgp2] - finf )
)/pow(r[xgp1],npow);
}
}
}
/* Upper x-bound: x(nx,:,:) --> xgp[xgp0] */
if (doBC[1]==1){
for (k=0;k<lsh[2];k++) {
for (j=0;j<lsh[1];j++) {
xgp1 = CCTK_GFINDEX3D(GH,lsh[0]-1,j,k);
xgp2 = CCTK_GFINDEX3D(GH,lsh[0]-2,j,k);
xgp3 = CCTK_GFINDEX3D(GH,lsh[0]-3,j,k);
xgp4 = CCTK_GFINDEX3D(GH,lsh[0]-4,j,k);
var[xgp1] = finf + (
pow(r[xgp3],npow)*( var[xgp4] - finf )
-3.0*pow(r[xgp3],npow)*( var[xgp3] - finf )
+3.0*pow(r[xgp2],npow)*( var[xgp2] - finf )
)/pow(r[xgp1],npow);
}
}
}
/* Lower y-bound */
if (doBC[2] == 1) {
for (k=0;k<lsh[2];k++) {
for (i=0;i<lsh[0];i++) {
ygp1 = CCTK_GFINDEX3D(GH,i,0,k);
ygp2 = CCTK_GFINDEX3D(GH,i,1,k);
ygp3 = CCTK_GFINDEX3D(GH,i,2,k);
ygp4 = CCTK_GFINDEX3D(GH,i,3,k);
var[ygp1] = finf + (
pow(r[ygp4],npow)*( var[ygp4] - finf )
-3.0*pow(r[ygp3],npow)*( var[ygp3] - finf )
+3.0*pow(r[ygp2],npow)*( var[ygp2] - finf )
)/pow(r[ygp1],npow);
}
}
}
/* Upper y bound */
if (doBC[3] == 1) {
for (k=0;k<lsh[2];k++) {
for (i=0;i<lsh[0];i++) {
ygp1 = CCTK_GFINDEX3D(GH,i,lsh[1]-1,k);
ygp2 = CCTK_GFINDEX3D(GH,i,lsh[1]-2,k);
ygp3 = CCTK_GFINDEX3D(GH,i,lsh[1]-3,k);
ygp4 = CCTK_GFINDEX3D(GH,i,lsh[1]-4,k);
var[ygp1] = finf + (
pow(r[ygp3],npow)*( var[ygp4] - finf )
-3.0*pow(r[ygp3],npow)*( var[ygp3] - finf )
+3.0*pow(r[ygp2],npow)*( var[ygp2] - finf )
)/pow(r[ygp1],npow);
}
}
}
/* Lower z-bound */
if (doBC[4]==1) {
for (j=0;j<lsh[1];j++) {
for (i=0;i<lsh[0];i++) {
zgp1 = CCTK_GFINDEX3D(GH,i,j,1);
zgp2 = CCTK_GFINDEX3D(GH,i,j,2);
zgp3 = CCTK_GFINDEX3D(GH,i,j,3);
zgp4 = CCTK_GFINDEX3D(GH,i,j,4);
var[zgp1] = finf + (
pow(r[zgp4],npow)*( var[zgp4] - finf )
-3.0*pow(r[zgp3],npow)*( var[zgp3] - finf )
+3.0*pow(r[zgp2],npow)*( var[zgp2] - finf )
)/pow(r[zgp1],npow);
}
}
}
/* Upper z-bound */
if (doBC[5] == 1) {
for (j=0;j<lsh[1];j++) {
for (i=0;i<lsh[0];i++) {
zgp1 = CCTK_GFINDEX3D(GH,i,j,lsh[2]-1);
zgp2 = CCTK_GFINDEX3D(GH,i,j,lsh[2]-2);
zgp3 = CCTK_GFINDEX3D(GH,i,j,lsh[2]-3);
zgp4 = CCTK_GFINDEX3D(GH,i,j,lsh[2]-4);
var[zgp1] = finf + (
pow(r[zgp3],npow)*( var[zgp4] - finf )
-3.0*pow(r[zgp3],npow)*( var[zgp3] - finf )
+3.0*pow(r[zgp2],npow)*( var[zgp2] - finf )
)/pow(r[zgp1],npow);
}
}
}
}
|
