blob: dcc63adda44e35b344d8126b8488cd205a05d40c (
plain)
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
|
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#if 1
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#endif
#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b))
static inline void steep(double *x, double *dx, double* dmx, const int i) {
if ( (x[i+1] - x[i]) * (x[i]-x[i-1]) > 0.0 ) {
dmx[i] = copysign(1.0,dx[i]) * MIN(fabs(dx[i]),
MIN(2.0*fabs(x[i]-x[i-1]),
2.0*fabs(x[i+1]-x[i])));
} else {
dmx[i] = 0.0;
}
}
static inline double approx_at_cell_interface(double* a, const int i) {
return 7.0/12.0*(a[i]+a[i+1]) - 1.0/12.0*(a[i-1]+a[i+2]);
}
static inline void monotonize(double* restrict xminus,
double* restrict x,
double* restrict xplus,
const int i) {
if ( !(xplus[i]==x[i] && x[i]==xminus[i])
&& ( (xplus[i]-x[i])*(x[i]-xminus[i]) <= 0.0 ) )
{
xminus[i] = x[i];
xplus[i] = x[i];
} else if( 6.0 * (xplus[i]-xminus[i]) *
(x[i]-0.5*(xplus[i]+xminus[i])) >
(xplus[i]-xminus[i])*(xplus[i]-xminus[i]) )
{
xminus[i] = 3.0*x[i]-2.0*xplus[i];
} else if( 6.0 * (xplus[i]-xminus[i]) *
(x[i]-0.5*(xplus[i]+xminus[i])) <
-(xplus[i]-xminus[i])*(xplus[i]-xminus[i]) )
{
xplus[i] = 3.0*x[i]-2.0*xminus[i];
}
return;
}
|
