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
|
#include <math.h>
#include <stddef.h>
#include <stdint.h>
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "Symmetry.h"
#define SQR(x) ((x) * (x))
void ns_mol_register(CCTK_ARGUMENTS)
{
MoLRegisterEvolved(CCTK_VarIndex("NullSurf::F"), CCTK_VarIndex("NullSurf::F_rhs"));
}
void ns_initial(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
for (int j = 0; j < cctkGH->cctk_lsh[2]; j++) {
for (int i = 0; i < cctkGH->cctk_lsh[0]; i++) {
const int idx = CCTK_GFINDEX3D(cctkGH, i, 0, j);
const double r2 = SQR(x[idx]) + SQR(z[idx]);
F[idx] = r2;
}
}
}
#define FD4(arr, stride, h_inv) \
((-1.0 * arr[2 * stride] + 8.0 * arr[1 * stride] - 8.0 * arr[-1 * stride] + 1.0 * arr[-2 * stride]) * h_inv / 12.0)
#define FD24(arr, stride, h_inv) \
((-1.0 * arr[2 * stride] + 16.0 * arr[1 * stride] - 30.0 * arr[0] + 16.0 * arr[-1 * stride] - 1.0 * arr[-2 * stride]) * SQR(h_inv) / 12.0)
static double diff4_dx(const double *arr, ptrdiff_t stride, double h_inv)
{
return FD4(arr, 1, h_inv);
}
static double diff4_dz(const double *arr, ptrdiff_t stride, double h_inv)
{
return FD4(arr, stride, h_inv);
}
void ns_eval_rhs(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
const int ghosts_x = cctkGH->cctk_nghostzones[0];
const int ghosts_z = cctkGH->cctk_nghostzones[2];
const ptrdiff_t stride_z = CCTK_GFINDEX3D(cctkGH, 0, 0, 1);
const double dx = x[1] - x[0];
const double dz = z[stride_z] - z[0];
const double dx_inv = 1.0 / dx;
const double dz_inv = 1.0 / dz;
for (int idx_z = ghosts_z; idx_z < cctkGH->cctk_lsh[2] - ghosts_z; idx_z++) {
for (int idx_x = ghosts_x; idx_x < cctkGH->cctk_lsh[0] - ghosts_x; idx_x++) {
const int idx = CCTK_GFINDEX3D(cctkGH, idx_x, 0, idx_z);
const double dF[3] = {diff4_dx(&F[idx], 1, dx_inv), 0.0, diff4_dz(&F[idx], stride_z, dz_inv)};
const double alpha_val = alpha[idx];
const double beta_val[3] = { beta1[idx], 0.0, beta3[idx] };
const double phi_val = phi[idx];
const double gtxx = gt11[idx];
const double gtxy = gt12[idx];
const double gtxz = gt13[idx];
const double gtyy = gt22[idx];
const double gtyz = gt23[idx];
const double gtzz = gt33[idx];
const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
double gtu[3][3], g4u[4][4];
double shift_term;
// \tilde{γ}^{ij}
gtu[0][0] = (gtyy * gtzz - SQR(gtyz)) / det;
gtu[1][1] = (gtxx * gtzz - SQR(gtxz)) / det;
gtu[2][2] = (gtxx * gtyy - SQR(gtxy)) / det;
gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det;
gtu[0][2] = (gtxy * gtyz - gtyy * gtxz) / det;
gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det;
gtu[1][0] = gtu[0][1];
gtu[2][0] = gtu[0][2];
gtu[2][1] = gtu[1][2];
g4u[0][0] = - 1.0 / SQR(alpha_val);
for (int i = 0; i < 3; i++) {
const double val = beta_val[i] / SQR(alpha_val);
g4u[1 + i][0] = val;
g4u[0][1 + i] = val;
}
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) {
g4u[1 + i][1 + j] = SQR(phi_val) * gtu[i][j] - beta_val[i] * beta_val[j] / SQR(alpha_val);
}
shift_term = 0.0;
for (int i = 0; i < 3; i++)
shift_term += g4u[0][1 + i] * dF[i];
double dF2 = 0.0;
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
dF2 += g4u[1 + i][1 + j] * dF[i] * dF[j];
F_rhs[idx] = (-shift_term + sqrt(SQR(shift_term) - g4u[0][0] * dF2)) / g4u[0][0];
}
}
}
void ns_register_symmetries(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
int sym[3] = { 1, 1, 1 };
int ret;
ret = SetCartSymVN(cctkGH, sym, "NullSurf::F");
if (ret != 0)
CCTK_WARN(0, "Error registering symmetries");
}
void ns_select_bc(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
int ret = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, "NullSurf::F", "none");
if (ret != 0)
CCTK_WARN(0, "Error registering boundaries");
}
|
