inoisyB/model_uniform.c at master · AFD-Illinois/inoisyB · GitHub

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

/*

noisy:

Generate a gaussian random field in 2D with
locally anisotropic correlation function,
locally varying correlation time.

Follows the technique of
Lindgren, Rue, and Lindstr\:om 2011, J.R. Statist. Soc. B 73, pp 423-498.
https://rss.onlinelibrary.wiley.com/doi/10.1111/j.1467-9868.2011.00777.x
in particular, implements eq. 17, which has power spectrum given by eq. 18.

Based on work by
David Daeyoung Lee
Charles Gammie
on applications in disk turbulence.

CFG 22 Dec 2019

*/

#include "noisy.h"

/*** this next bit describes the model ***/

double correlation_length(double x,double y)
{

    return( 1.0 );

}

double correlation_time(double x,double y)
{

    return( 1.0 );
}

double diffusion_coefficient(double x,double y)
{

    double K = 1.0;

    return( K );
}

/* return principal axes of diffusion tensor */

void principal_axis_func(double x, double y, double *e1x, double *e1y, double *e2x, double *e2y)
{
    double phi,s,c;

    phi = M_PI/4.;

    c = cos(phi);
    s = sin(phi);

    *e1x = c;
    *e1y = s;
    *e2x = -s;
    *e2y = c;
}

/* return advection velocity as a function of position */

void advection_velocity(double x, double y, double va[2])
{

    double vx = 0. ;
    double vy = 0. ;

    va[0] = vx ;
    va[1] = vy ;
}


/* return envelope function; image is exp(-del)*envelope */

double envelope(double x, double y)
{
    return( 1.0 );
}

/* noise model */
void noise_model(double del_noise[][N], double dt)
{

    int i,j;
    static int first_call = 1;
    static gsl_rng * r;

#if 0
	/* white noise in time model */
    if(first_call) {
        r = gsl_rng_alloc(gsl_rng_mt19937); /* Mersenne twister */
        gsl_rng_set(r, 0);
        first_call = 0;
    }

    /* update del */
    for(i=0;i<N;i++)
    for(j=0;j<N;j++) {
        del[i][j] += PARAM_EPS*gsl_ran_gaussian_ziggurat(r,1.0);
    }
#endif

#if 1
	/* static model */
	static double save_noise[N][N];
    if(first_call) {
        r = gsl_rng_alloc(gsl_rng_mt19937); /* Mersenne twister */
        gsl_rng_set(r, 0);

		for(i=0;i<N;i++)
		for(j=0;j<N;j++) {
			save_noise[i][j] = PARAM_EPS*gsl_ran_gaussian_ziggurat(r,1.0);
		}

        first_call = 0;
	}
    for(i=0;i<N;i++)
    for(j=0;j<N;j++) {
		del_noise[i][j] = dt*save_noise[i][j];
	}
#endif

}