2dnavierstokes/init_flow.c at main · martinjamesphd/2dnavierstokes · 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
#include<math.h>
#include<stdio.h>
#include<rfftw_threads.h>
#include<time.h>
#include"ns2d.h"

//function to initialize time increment and density of states
void init_den_state(int *den_state)
{
	int i, j, k, ij;
	int y_size_f=sys_size/2+1;
	int y_alias = sys_size/3;

	for(i=0;i<y_size_f;i++)
		den_state[i]=0;

	//find density of state for each k
	for(i=0;i<y_size_f;i++)
		for(j=0;j<y_size_f;j++)
		{
			k=(int)sqrt(i*i+j*j);
			if(k<y_size_f) den_state[k]++;
		}

	for(k=0;k<y_size_f;k++)
		den_state[k]*=4;	//accounts for each quadrant

	for(i=0;i<y_size_f;i++)
		for(j=0;j<y_size_f;j++)
		{
			ij=i*y_size_f+j;
			exp_dt[ij]=exp(-(vis*(i*i+j*j)+mu)*dt/2);

			if((i*i+j*j)>y_alias*y_alias)	//accounts for aliasing (vorticity evolution uses this factor)
				exp_dt[ij]=0.0;
		}

	for(i=y_size_f;i<sys_size;i++)	//i greater than N/2 corresponds to kx = i-N (where N is sys_size)
		for(j=0;j<y_size_f;j++)
		{
		 	ij=i*y_size_f+j;
		 	exp_dt[ij]=exp(-(vis*((sys_size-i)*(sys_size-i)+j*j)+mu)*dt/2);

		 	if(((sys_size-i)*(sys_size-i)+j*j)>y_alias*y_alias)
		 		exp_dt[ij]=0.0;
		}
}

//initializes omega in fourier space
void init_omega(fftw_complex *omega_ft, int *den_state)
{
	int i, j, k, ij, kSqr;
	int y_size_f = sys_size/2+1;

	double tempomega;

	srand(time(NULL));	//seeds random function using current time

	for(i=0;i<y_size_f;i++)
		for(j=0;j<y_size_f;j++)
		{
			ij=i*y_size_f+j;
			kSqr=i*i+j*j;
			k=(int)sqrt(kSqr);

			if(k<10&&k<y_size_f)
			{
				tempomega= kSqr*kSqr*exp(-kSqr*kSqr/2);
				tempomega/=den_state[k];
				tempomega*=(double)sys_size*sys_size;

				omega_ft[ij].re=tempomega*cos(2.0*M_PI*rand()/RAND_MAX);
				omega_ft[ij].im=tempomega*sin(2.0*M_PI*rand()/RAND_MAX);
			}
			else
			{
				omega_ft[ij].re=0;
				omega_ft[ij].im=0;
			}
		}
	for(i=y_size_f;i<sys_size;i++)//i greater than N/2 corresponds to kx = i-N (where N is sys_size)
		for(j=0;j<y_size_f;j++)
	 	{
	 		ij=i*y_size_f+j;
			kSqr=(sys_size-i)*(sys_size-i)+j*j;
			k=(int)sqrt(kSqr);

			if(k<10)
			{
				tempomega= kSqr*kSqr*exp(-kSqr*kSqr/2);
				tempomega/=den_state[k];
				tempomega*=(double)sys_size*sys_size;

				omega_ft[ij].re=tempomega*cos(2.0*M_PI*rand()/RAND_MAX);
				omega_ft[ij].im=tempomega*sin(2.0*M_PI*rand()/RAND_MAX);
			}
			else
			{
				omega_ft[ij].re=0;
				omega_ft[ij].im=0;
			}
		}
}

//Generates stirring fomega
void gen_force(double famp, int kf)
{
	int y_size_r=sys_size+2;
	int i, j, ij;

	double *fomega = (double*)fomega_ft;

	for(i=0;i<sys_size;i++)
		for(j=0;j<y_size_r;j++)
		{
			ij = i*y_size_r+j;
			fomega[ij]=-kf*famp*cos(kf*2.0*M_PI*i/(sys_size));
		}
	rfftwnd_threads_one_real_to_complex(nthreads, p_for, fomega, 0);

}