forked from AFD-Illinois/igrmonty2d
-
Notifications
You must be signed in to change notification settings - Fork 0
/
init_geometry.c
66 lines (43 loc) · 1.39 KB
/
init_geometry.c
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
#include "decs.h"
/*
set up the metric (indicies both up and down), and the
connection coefficients on a grid, based
on a HARM simulation grid.
In principle the geometry could be sampled on a finer or
coarser grid than the simulation grid.
These routines are taken directly out of HARM.
They require the code be compiled against the
Gnu scientific library (GSL).
CFG 21 July 06
*/
gsl_matrix *gsl_gcov, *gsl_gcon;
gsl_permutation *perm;
#pragma omp threadprivate (gsl_gcov, gsl_gcon, perm)
/* assumes gcov has been set first; returns determinant */
double gdet_func(double gcov[][NDIM])
{
double d;
int signum;
if (gsl_gcov == NULL) {
gsl_gcov = gsl_matrix_alloc(NDIM, NDIM);
gsl_gcon = gsl_matrix_alloc(NDIM, NDIM);
perm = gsl_permutation_alloc(NDIM);
}
DLOOP gsl_matrix_set(gsl_gcov, k, l, gcov[k][l]);
gsl_linalg_LU_decomp(gsl_gcov, perm, &signum);
d = gsl_linalg_LU_det(gsl_gcov, signum);
return (sqrt(fabs(d)));
}
void gcon_func(double gcov[][NDIM], double gcon[][NDIM])
{
int signum;
if (gsl_gcov == NULL) {
gsl_gcov = gsl_matrix_alloc(NDIM, NDIM);
gsl_gcon = gsl_matrix_alloc(NDIM, NDIM);
perm = gsl_permutation_alloc(NDIM);
}
DLOOP gsl_matrix_set(gsl_gcov, k, l, gcov[k][l]);
gsl_linalg_LU_decomp(gsl_gcov, perm, &signum);
gsl_linalg_LU_invert(gsl_gcov, perm, gsl_gcon);
DLOOP gcon[k][l] = gsl_matrix_get(gsl_gcon, k, l);
}