Initial commit
This commit is contained in:
Ian Jauslin
282
src/main.c
Normal file
282
src/main.c
Normal file
@ -0,0 +1,282 @@
|
||||
#define VERSION "0.0"
|
||||
|
||||
#include <math.h>
|
||||
#include <complex.h>
|
||||
#include <fftw3.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include "navier-stokes.h"
|
||||
|
||||
// usage message
|
||||
int print_usage();
|
||||
// read command line arguments
|
||||
int read_args(int argc, const char* argv[], ns_params* params, unsigned int* nsteps, unsigned int* computation_nr);
|
||||
|
||||
// compute enstrophy as a function of time in the I-NS equation
|
||||
int enstrophy(ns_params params, unsigned int Nsteps);
|
||||
|
||||
|
||||
#define COMPUTATION_ENSTROPHY 1
|
||||
int main (int argc, const char* argv[]){
|
||||
ns_params params;
|
||||
unsigned int nsteps;
|
||||
int ret;
|
||||
unsigned int computation_nr;
|
||||
|
||||
// default computation: phase diagram
|
||||
computation_nr=COMPUTATION_ENSTROPHY;
|
||||
|
||||
// read command line arguments
|
||||
ret=read_args(argc, argv, ¶ms, &nsteps, &computation_nr);
|
||||
if(ret<0){
|
||||
return(-1);
|
||||
}
|
||||
if(ret>0){
|
||||
return(0);
|
||||
}
|
||||
|
||||
// enstrophy
|
||||
if(computation_nr==COMPUTATION_ENSTROPHY){
|
||||
enstrophy(params, nsteps);
|
||||
}
|
||||
|
||||
return(0);
|
||||
}
|
||||
|
||||
// usage message
|
||||
int print_usage(){
|
||||
fprintf(stderr, "usage:\n nstrophy enstrophy [-h timestep] [-K modes] [-v] [-N nsteps]\n\n nstrophy -V [-v]\n\n");
|
||||
return(0);
|
||||
}
|
||||
|
||||
// read command line arguments
|
||||
#define CP_FLAG_TIMESTEP 1
|
||||
#define CP_FLAG_NSTEPS 2
|
||||
#define CP_FLAG_MODES 2
|
||||
#define CP_FLAG_NU 3
|
||||
int read_args(int argc, const char* argv[], ns_params* params, unsigned int* nsteps, unsigned int* computation_nr){
|
||||
int i;
|
||||
int ret;
|
||||
// temporary int
|
||||
int tmp_int;
|
||||
// temporary unsigned int
|
||||
unsigned int tmp_uint;
|
||||
// temporary double
|
||||
double tmp_double;
|
||||
// pointers
|
||||
char* ptr;
|
||||
// flag that indicates what argument is being read
|
||||
int flag=0;
|
||||
// print version and exit
|
||||
char Vflag=0;
|
||||
|
||||
// defaults
|
||||
/*
|
||||
params->h=6.103515625e-05;
|
||||
params->K=16;
|
||||
*nsteps=16777216;
|
||||
params->nu=4.9632717887631524e-05;
|
||||
*/
|
||||
params->h=1e-5;
|
||||
params->K=16;
|
||||
*nsteps=10000000;
|
||||
params->nu=1e-4;
|
||||
|
||||
// loop over arguments
|
||||
for(i=1;i<argc;i++){
|
||||
// flag
|
||||
if(argv[i][0]=='-'){
|
||||
for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
|
||||
switch(*ptr){
|
||||
// timestep
|
||||
case 'h':
|
||||
flag=CP_FLAG_TIMESTEP;
|
||||
break;
|
||||
// nsteps
|
||||
case 'N':
|
||||
flag=CP_FLAG_NSTEPS;
|
||||
break;
|
||||
// modes
|
||||
case 'K':
|
||||
flag=CP_FLAG_MODES;
|
||||
break;
|
||||
// friction
|
||||
case 'n':
|
||||
flag=CP_FLAG_NU;
|
||||
break;
|
||||
// print version
|
||||
case 'V':
|
||||
Vflag=1;
|
||||
break;
|
||||
default:
|
||||
fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
|
||||
print_usage();
|
||||
return(-1);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
// timestep
|
||||
else if(flag==CP_FLAG_TIMESTEP){
|
||||
ret=sscanf(argv[i],"%lf",&tmp_double);
|
||||
if(ret!=1){
|
||||
fprintf(stderr, "error: '-h' should be followed by a double\n got '%s'\n",argv[i]);
|
||||
return(-1);
|
||||
}
|
||||
params->h=tmp_double;
|
||||
flag=0;
|
||||
}
|
||||
// nsteps
|
||||
else if(flag==CP_FLAG_NSTEPS){
|
||||
ret=sscanf(argv[i],"%u",&tmp_uint);
|
||||
if(ret!=1){
|
||||
fprintf(stderr, "error: '-N' should be followed by an unsigned int\n got '%s'\n",argv[i]);
|
||||
return(-1);
|
||||
}
|
||||
*nsteps=tmp_uint;
|
||||
flag=0;
|
||||
}
|
||||
// nsteps
|
||||
else if(flag==CP_FLAG_MODES){
|
||||
ret=sscanf(argv[i],"%d",&tmp_int);
|
||||
if(ret!=1){
|
||||
fprintf(stderr, "error: '-K' should be followed by an int\n got '%s'\n",argv[i]);
|
||||
return(-1);
|
||||
}
|
||||
params->K=tmp_uint;
|
||||
flag=0;
|
||||
}
|
||||
// friction
|
||||
else if(flag==CP_FLAG_TIMESTEP){
|
||||
ret=sscanf(argv[i],"%lf",&tmp_double);
|
||||
if(ret!=1){
|
||||
fprintf(stderr, "error: '-n' should be followed by a double\n got '%s'\n",argv[i]);
|
||||
return(-1);
|
||||
}
|
||||
params->nu=tmp_double;
|
||||
flag=0;
|
||||
}
|
||||
// computation to run
|
||||
else{
|
||||
if(strcmp(argv[i], "enstrophy")==0){
|
||||
*computation_nr=COMPUTATION_ENSTROPHY;
|
||||
}
|
||||
else{
|
||||
fprintf(stderr, "error: unrecognized computation: '%s'\n",argv[i]);
|
||||
print_usage();
|
||||
return(-1);
|
||||
}
|
||||
flag=0;
|
||||
}
|
||||
}
|
||||
|
||||
// print version and exit
|
||||
if(Vflag==1){
|
||||
printf("nstrophy " VERSION "\n");
|
||||
return(1);
|
||||
}
|
||||
|
||||
return(0);
|
||||
}
|
||||
|
||||
// compute enstrophy as a function of time in the I-NS equation
|
||||
int enstrophy(ns_params params, unsigned int Nsteps){
|
||||
_Complex double* u;
|
||||
_Complex double* tmp1;
|
||||
_Complex double* tmp2;
|
||||
_Complex double* tmp3;
|
||||
_Complex double alpha;
|
||||
_Complex double avg;
|
||||
unsigned int t;
|
||||
int kx,ky;
|
||||
fft_vects fft_vects;
|
||||
|
||||
// sizes
|
||||
params.S=2*params.K+1;
|
||||
params.N=4*params.K+1;
|
||||
|
||||
// velocity field
|
||||
u=calloc(sizeof(_Complex double),params.S*params.S);
|
||||
params.g=calloc(sizeof(_Complex double),params.S*params.S);
|
||||
// allocate tmp vectors for computation
|
||||
tmp1=calloc(sizeof(_Complex double),params.S*params.S);
|
||||
tmp2=calloc(sizeof(_Complex double),params.S*params.S);
|
||||
tmp3=calloc(sizeof(_Complex double),params.S*params.S);
|
||||
|
||||
// initial value
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
//u[KLOOKUP(kx,ky,params.S)]=kx*kx*ky*ky*exp(-(kx*kx+ky*ky));
|
||||
if((kx==1 && ky==0) || (kx==-1 && ky==0)){
|
||||
u[KLOOKUP(kx,ky,params.S)]=1;
|
||||
}
|
||||
else{
|
||||
u[KLOOKUP(kx,ky,params.S)]=0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// driving force
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
//params.g[KLOOKUP(kx,ky,params.S)]=sqrt(kx*kx*ky*ky)*exp(-(kx*kx+ky*ky));
|
||||
if((kx==2 && ky==-1) || (kx==-2 && ky==1)){
|
||||
params.g[KLOOKUP(kx,ky,params.S)]=1;
|
||||
}
|
||||
else{
|
||||
params.g[KLOOKUP(kx,ky,params.S)]=0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// prepare vectors for fft
|
||||
fft_vects.fft1=fftw_malloc(sizeof(fftw_complex)*params.N*params.N);
|
||||
fft_vects.fft1_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.fft1, fft_vects.fft1, FFTW_FORWARD, FFTW_MEASURE);
|
||||
fft_vects.fft2=fftw_malloc(sizeof(fftw_complex)*params.N*params.N);
|
||||
fft_vects.fft2_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.fft2, fft_vects.fft2, FFTW_FORWARD, FFTW_MEASURE);
|
||||
fft_vects.invfft=fftw_malloc(sizeof(fftw_complex)*params.N*params.N);
|
||||
fft_vects.invfft_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.invfft, fft_vects.invfft, FFTW_BACKWARD, FFTW_MEASURE);
|
||||
|
||||
// init running average
|
||||
avg=0;
|
||||
|
||||
// iterate
|
||||
for(t=0;t<Nsteps;t++){
|
||||
ins_step(u, params, fft_vects, tmp1, tmp2, tmp3);
|
||||
alpha=compute_alpha(u, params);
|
||||
|
||||
// to avoid errors building up in imaginary part
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
u[KLOOKUP(kx,ky,params.S)]=__real__ u[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
|
||||
// running average
|
||||
if(t>0){
|
||||
avg=avg-(avg-alpha)/t;
|
||||
}
|
||||
|
||||
if(t>0 && t%1000==0){
|
||||
fprintf(stderr,"%8d % .8e % .8e % .8e % .8e\n",t, __real__ avg, __imag__ avg, __real__ alpha, __imag__ alpha);
|
||||
printf("%8d % .8e % .8e % .8e % .8e\n",t, __real__ avg, __imag__ avg, __real__ alpha, __imag__ alpha);
|
||||
}
|
||||
}
|
||||
|
||||
// free memory
|
||||
fftw_destroy_plan(fft_vects.fft1_plan);
|
||||
fftw_destroy_plan(fft_vects.fft2_plan);
|
||||
fftw_destroy_plan(fft_vects.invfft_plan);
|
||||
fftw_free(fft_vects.fft1);
|
||||
fftw_free(fft_vects.fft2);
|
||||
fftw_free(fft_vects.invfft);
|
||||
|
||||
free(tmp3);
|
||||
free(tmp2);
|
||||
free(tmp1);
|
||||
free(params.g);
|
||||
free(u);
|
||||
|
||||
return(0);
|
||||
}
|
||||
299
src/navier-stokes.c
Normal file
299
src/navier-stokes.c
Normal file
@ -0,0 +1,299 @@
|
||||
#include "navier-stokes.h"
|
||||
#include <math.h>
|
||||
|
||||
#define M_PI 3.14159265358979323846
|
||||
|
||||
#define CHK 1
|
||||
|
||||
// next time step for Irreversible Navier-Stokes equation
|
||||
int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3){
|
||||
int kx,ky;
|
||||
|
||||
// k1
|
||||
ins_rhs(tmp1, u, params, vects);
|
||||
// add to output
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
tmp3[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/6*tmp1[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
|
||||
// u+h*k1/2
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/2*tmp1[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
// k2
|
||||
ins_rhs(tmp1, tmp2, params, vects);
|
||||
// add to output
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
tmp3[KLOOKUP(kx,ky,params.S)]+=params.h/3*tmp1[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
|
||||
// u+h*k2/2
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/2*tmp1[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
// k3
|
||||
ins_rhs(tmp1, tmp2, params, vects);
|
||||
// add to output
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
tmp3[KLOOKUP(kx,ky,params.S)]+=params.h/3*tmp1[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
|
||||
// u+h*k3
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h*tmp1[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
// k4
|
||||
ins_rhs(tmp1, tmp2, params, vects);
|
||||
// add to output
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
u[KLOOKUP(kx,ky,params.S)]=tmp3[KLOOKUP(kx,ky,params.S)]+params.h/6*tmp1[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
|
||||
return(0);
|
||||
}
|
||||
|
||||
// right side of Irreversible Navier-Stokes equation
|
||||
int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vects vects){
|
||||
int kx,ky;
|
||||
|
||||
#if CHK==0
|
||||
// F(u/|p|)*F(q1*q2*u/|q|)
|
||||
// init to 0
|
||||
for(kx=0; kx<params.N*params.N; kx++){
|
||||
vects.fft1[kx]=0;
|
||||
vects.fft2[kx]=0;
|
||||
vects.invfft[kx]=0;
|
||||
}
|
||||
// fill modes
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
if(kx!=0 || ky!=0){
|
||||
vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
|
||||
vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*ky*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
|
||||
}
|
||||
}
|
||||
}
|
||||
// fft
|
||||
fftw_execute(vects.fft1_plan);
|
||||
fftw_execute(vects.fft2_plan);
|
||||
// write to invfft
|
||||
for(kx=-2*params.K;kx<=2*params.K;kx++){
|
||||
for(ky=-2*params.K;ky<=2*params.K;ky++){
|
||||
vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
|
||||
}
|
||||
}
|
||||
// inverse fft
|
||||
fftw_execute(vects.invfft_plan);
|
||||
|
||||
// write out
|
||||
for(kx=0; kx<params.S*params.S; kx++){
|
||||
out[kx]=0;
|
||||
}
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
if(kx!=0 || ky!=0){
|
||||
out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*kx-ky*ky)/params.N/params.N;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// F(u/|p|)*F((q1*q1-q2*q2)*u/|q|)
|
||||
// init to 0
|
||||
for(kx=0; kx<params.N*params.N; kx++){
|
||||
vects.fft2[kx]=0;
|
||||
vects.invfft[kx]=0;
|
||||
}
|
||||
// fill modes
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
if(kx!=0 || ky!=0){
|
||||
vects.fft2[KLOOKUP(kx,ky,params.N)]=(kx*kx-ky*ky)*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
|
||||
}
|
||||
}
|
||||
}
|
||||
// fft
|
||||
fftw_execute(vects.fft2_plan);
|
||||
// write to invfft
|
||||
for(kx=-2*params.K;kx<=2*params.K;kx++){
|
||||
for(ky=-2*params.K;ky<=2*params.K;ky++){
|
||||
vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
|
||||
}
|
||||
}
|
||||
// inverse fft
|
||||
fftw_execute(vects.invfft_plan);
|
||||
|
||||
// write out
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
if(kx!=0 || ky!=0){
|
||||
out[KLOOKUP(kx,ky,params.S)]=out[KLOOKUP(kx,ky,params.S)]-4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*ky)/params.N/params.N;
|
||||
}
|
||||
}
|
||||
}
|
||||
#elif CHK == 1
|
||||
// F(-p2/|p|*u)*F(q1*|q|*u)
|
||||
// init to 0
|
||||
for(kx=0; kx<params.N*params.N; kx++){
|
||||
vects.fft1[kx]=0;
|
||||
vects.fft2[kx]=0;
|
||||
vects.invfft[kx]=0;
|
||||
}
|
||||
// fill modes
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
if(kx!=0 || ky!=0){
|
||||
vects.fft1[KLOOKUP(kx,ky,params.N)]=-kx/sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
|
||||
vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
}
|
||||
// fft
|
||||
fftw_execute(vects.fft1_plan);
|
||||
fftw_execute(vects.fft2_plan);
|
||||
// write to invfft
|
||||
for(kx=-2*params.K;kx<=2*params.K;kx++){
|
||||
for(ky=-2*params.K;ky<=2*params.K;ky++){
|
||||
vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)] - vects.fft1[KLOOKUP(ky,kx,params.N)]*vects.fft2[KLOOKUP(ky,kx,params.N)];
|
||||
}
|
||||
}
|
||||
|
||||
// inverse fft
|
||||
fftw_execute(vects.invfft_plan);
|
||||
|
||||
|
||||
// write out
|
||||
for(kx=0; kx<params.S*params.S; kx++){
|
||||
out[kx]=0;
|
||||
}
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
if(kx!=0 || ky!=0){
|
||||
out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu/params.S*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+2*M_PI/sqrt(kx*kx+ky*ky)/params.S*vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#elif CHK==2
|
||||
// F(u)*F(q1*u)
|
||||
// init to 0
|
||||
for(kx=0; kx<params.N*params.N; kx++){
|
||||
vects.fft1[kx]=0;
|
||||
vects.fft2[kx]=0;
|
||||
vects.invfft[kx]=0;
|
||||
}
|
||||
// fill modes
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
// u_k
|
||||
vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)];
|
||||
// 2i\pi k_x u_k
|
||||
vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*kx*u[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
// fft
|
||||
fftw_execute(vects.fft1_plan);
|
||||
fftw_execute(vects.fft2_plan);
|
||||
// write to invfft
|
||||
for(kx=-2*params.K;kx<=2*params.K;kx++){
|
||||
for(ky=-2*params.K;ky<=2*params.K;ky++){
|
||||
vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
|
||||
}
|
||||
}
|
||||
|
||||
// F(p1/p2*u)*F(q2*u)
|
||||
// init to 0
|
||||
for(kx=0; kx<params.N*params.N; kx++){
|
||||
vects.fft1[kx]=0;
|
||||
vects.fft2[kx]=0;
|
||||
}
|
||||
// fill modes
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
// k_x/k_y u_k
|
||||
if(ky!=0){
|
||||
vects.fft1[KLOOKUP(kx,ky,params.N)]=kx/ky*u[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
// 2i\pi k_y u_k
|
||||
vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*ky*u[KLOOKUP(kx,ky,params.S)];
|
||||
}
|
||||
}
|
||||
// fft
|
||||
fftw_execute(vects.fft1_plan);
|
||||
fftw_execute(vects.fft2_plan);
|
||||
// write to invfft
|
||||
for(kx=-2*params.K;kx<=2*params.K;kx++){
|
||||
for(ky=-2*params.K;ky<=2*params.K;ky++){
|
||||
vects.invfft[KLOOKUP(kx,ky,params.N)]+=-vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
|
||||
}
|
||||
}
|
||||
|
||||
// inverse fft
|
||||
fftw_execute(vects.invfft_plan);
|
||||
|
||||
/*
|
||||
// check: convolution instead of fft
|
||||
for(kx=0; kx<params.S*params.S; kx++){
|
||||
out[kx]=0;
|
||||
}
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
for(px=-params.K;px<=params.K;px++){
|
||||
for(py=-params.K;py<=params.K;py++){
|
||||
if(kx-px<=params.K && kx-px>=-params.K && ky-py<=params.K && ky-py>=-params.K){
|
||||
out[KLOOKUP(kx,ky,params.S)]+=2*M_PI*I*(u[KLOOKUP(px,py,params.S)]*(kx-px)*u[KLOOKUP(kx-px,ky-py,params.S)]-(py==0?0:px/py*u[KLOOKUP(px,py,params.S)]*(ky-py)*u[KLOOKUP(kx-px,ky-py,params.S)]));
|
||||
}
|
||||
}
|
||||
}
|
||||
dd=(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])*(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])+(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)])*(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)]);
|
||||
if(dd>1e-25){
|
||||
printf("%d %d % .8e\n",kx,ky, dd);
|
||||
}
|
||||
}
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
// write out
|
||||
for(kx=0; kx<params.S*params.S; kx++){
|
||||
out[kx]=0;
|
||||
}
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
return(0);
|
||||
}
|
||||
|
||||
|
||||
// compute alpha
|
||||
_Complex double compute_alpha(_Complex double* u, ns_params params){
|
||||
_Complex double num=0;
|
||||
_Complex double denom=0;
|
||||
int kx,ky;
|
||||
|
||||
for(kx=-params.K;kx<=params.K;kx++){
|
||||
for(ky=-params.K;ky<=params.K;ky++){
|
||||
denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]*conj(u[KLOOKUP(kx,ky,params.S)])*(1+(ky!=0?kx*kx/ky/ky:0));
|
||||
num+=(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]*conj(params.g[KLOOKUP(kx,ky,params.S)])*(1+(ky!=0?kx*kx/ky/ky:0));
|
||||
}
|
||||
}
|
||||
|
||||
return(num/denom);
|
||||
}
|
||||
48
src/navier-stokes.h
Normal file
48
src/navier-stokes.h
Normal file
@ -0,0 +1,48 @@
|
||||
#ifndef NAVIERSTOKES_H
|
||||
#define NAVIERSTOKES_H
|
||||
|
||||
#include <complex.h>
|
||||
#include <fftw3.h>
|
||||
|
||||
// to extract elements from array of size S representing a function of momentum, use
|
||||
// array[KEXTRACT(kx,ky,size)]
|
||||
#define KLOOKUP(X,Y,S) (X>=0?X:S+X)*S+(Y>=0?Y:S+Y)
|
||||
|
||||
|
||||
// parameters for the NS equation
|
||||
typedef struct ns_params {
|
||||
// number of modes
|
||||
int K;
|
||||
// 2*K+1
|
||||
int S;
|
||||
// 4*K+1
|
||||
int N;
|
||||
// forcing term
|
||||
_Complex double* g;
|
||||
// time step
|
||||
double h;
|
||||
// friction
|
||||
double nu;
|
||||
} ns_params;
|
||||
|
||||
// arrays on which the ffts are performed
|
||||
typedef struct fft_vects {
|
||||
fftw_complex* fft1;
|
||||
fftw_complex* fft2;
|
||||
fftw_complex* invfft;
|
||||
fftw_plan fft1_plan;
|
||||
fftw_plan fft2_plan;
|
||||
fftw_plan invfft_plan;
|
||||
} fft_vects;
|
||||
|
||||
// next time step for Irreversible Navier-Stokes equation
|
||||
int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3);
|
||||
|
||||
// right side of Irreversible Navier-Stokes equation
|
||||
int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vects vects);
|
||||
|
||||
// compute alpha
|
||||
_Complex double compute_alpha(_Complex double* u, ns_params params);
|
||||
|
||||
#endif
|
||||
|
||||
Reference in New Issue
Block a user