2018-01-11 22:48:14 +00:00
|
|
|
#include "navier-stokes.h"
|
2022-05-27 16:09:17 -04:00
|
|
|
#include "io.h"
|
2018-01-11 22:48:14 +00:00
|
|
|
#include <math.h>
|
2022-05-18 09:57:06 +02:00
|
|
|
#include <stdlib.h>
|
2018-01-11 22:48:14 +00:00
|
|
|
|
2022-05-18 09:57:06 +02:00
|
|
|
// compute solution as a function of time
|
|
|
|
int uk(
|
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
int N1,
|
|
|
|
int N2,
|
|
|
|
unsigned int nsteps,
|
|
|
|
double nu,
|
|
|
|
double delta,
|
2022-05-25 11:12:02 -04:00
|
|
|
double L,
|
2022-05-26 15:05:30 -04:00
|
|
|
_Complex double* u0,
|
2022-05-26 15:16:44 -04:00
|
|
|
_Complex double* g,
|
2022-05-18 23:52:01 +02:00
|
|
|
unsigned int print_freq,
|
2022-05-27 17:06:38 -04:00
|
|
|
unsigned int starting_time,
|
2022-05-27 16:09:17 -04:00
|
|
|
unsigned int nthreads,
|
|
|
|
FILE* savefile
|
2022-05-18 09:57:06 +02:00
|
|
|
){
|
|
|
|
_Complex double* u;
|
|
|
|
_Complex double* tmp1;
|
|
|
|
_Complex double* tmp2;
|
|
|
|
_Complex double* tmp3;
|
|
|
|
unsigned int t;
|
|
|
|
fft_vect fft1;
|
|
|
|
fft_vect fft2;
|
|
|
|
fft_vect ifft;
|
|
|
|
int kx,ky;
|
|
|
|
|
2022-05-18 23:52:01 +02:00
|
|
|
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
|
2022-05-26 15:05:30 -04:00
|
|
|
// copy initial condition
|
|
|
|
copy_u(u, u0, K1, K2);
|
2022-05-18 09:57:06 +02:00
|
|
|
|
2022-05-19 17:51:45 +02:00
|
|
|
// print column headers
|
|
|
|
printf("# 1:i 2:t ");
|
|
|
|
t=3;
|
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for (ky=-K2;ky<=K2;ky++){
|
|
|
|
printf(" %6d:(%4d,%4d)r ",t,kx,ky);
|
|
|
|
t++;
|
|
|
|
printf(" %6d:(%4d,%4d)i ",t,kx,ky);
|
|
|
|
t++;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2022-05-18 09:57:06 +02:00
|
|
|
// iterate
|
2022-05-27 17:06:38 -04:00
|
|
|
for(t=starting_time;t<starting_time+nsteps;t++){
|
2022-05-25 11:12:02 -04:00
|
|
|
ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
|
2022-05-18 09:57:06 +02:00
|
|
|
|
|
|
|
if(t%print_freq==0){
|
2022-05-18 10:42:30 +02:00
|
|
|
fprintf(stderr,"%d % .8e ",t,t*delta);
|
|
|
|
printf("%8d % .15e ",t,t*delta);
|
2022-05-18 09:57:06 +02:00
|
|
|
|
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for (ky=-K2;ky<=K2;ky++){
|
|
|
|
if (kx*kx+ky*ky<=1){
|
|
|
|
fprintf(stderr,"% .8e % .8e ",__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)], __imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]);
|
|
|
|
|
|
|
|
}
|
2022-05-18 10:42:30 +02:00
|
|
|
printf("% .15e % .15e ",__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)], __imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]);
|
2022-05-18 09:57:06 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
fprintf(stderr,"\n");
|
|
|
|
printf("\n");
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2022-05-27 16:09:17 -04:00
|
|
|
// save final entry to savefile
|
|
|
|
write_u(u, K1, K2, savefile);
|
|
|
|
|
2022-05-18 09:57:06 +02:00
|
|
|
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
|
|
|
|
return(0);
|
|
|
|
}
|
|
|
|
|
2022-05-26 21:26:13 -04:00
|
|
|
// compute energy, enstrophy, alpha as a function of time in the I-NS equation
|
|
|
|
int eea(
|
2022-05-18 09:57:06 +02:00
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
int N1,
|
|
|
|
int N2,
|
|
|
|
unsigned int nsteps,
|
|
|
|
double nu,
|
|
|
|
double delta,
|
2022-05-25 11:12:02 -04:00
|
|
|
double L,
|
2022-05-26 15:05:30 -04:00
|
|
|
_Complex double* u0,
|
2022-05-26 15:16:44 -04:00
|
|
|
_Complex double* g,
|
2022-05-18 23:52:01 +02:00
|
|
|
unsigned int print_freq,
|
2022-05-27 17:06:38 -04:00
|
|
|
unsigned int starting_time,
|
2022-05-27 16:09:17 -04:00
|
|
|
unsigned int nthreads,
|
|
|
|
FILE* savefile
|
2022-05-18 09:57:06 +02:00
|
|
|
){
|
|
|
|
_Complex double* u;
|
|
|
|
_Complex double* tmp1;
|
|
|
|
_Complex double* tmp2;
|
|
|
|
_Complex double* tmp3;
|
2022-05-26 21:26:13 -04:00
|
|
|
double alpha, energy, enstrophy;
|
|
|
|
double avg_e,avg_a,avg_en;
|
2022-05-18 09:57:06 +02:00
|
|
|
unsigned int t;
|
|
|
|
fft_vect fft1;
|
|
|
|
fft_vect fft2;
|
|
|
|
fft_vect ifft;
|
|
|
|
|
2022-05-18 23:52:01 +02:00
|
|
|
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
|
2022-05-26 15:05:30 -04:00
|
|
|
// copy initial condition
|
|
|
|
copy_u(u, u0, K1, K2);
|
2022-05-18 09:57:06 +02:00
|
|
|
|
|
|
|
|
|
|
|
// init running average
|
2022-05-26 21:26:13 -04:00
|
|
|
avg_e=0;
|
|
|
|
avg_a=0;
|
|
|
|
avg_en=0;
|
2022-05-18 09:57:06 +02:00
|
|
|
|
|
|
|
// iterate
|
2022-05-27 17:06:38 -04:00
|
|
|
for(t=starting_time;t<starting_time+nsteps;t++){
|
2022-05-25 11:12:02 -04:00
|
|
|
ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
|
2022-05-26 21:26:13 -04:00
|
|
|
|
|
|
|
energy=compute_energy(u, K1, K2);
|
2022-05-18 09:57:06 +02:00
|
|
|
alpha=compute_alpha(u, K1, K2, g);
|
2022-05-26 21:26:13 -04:00
|
|
|
enstrophy=compute_enstrophy(u, K1, K2, L);
|
2022-05-18 09:57:06 +02:00
|
|
|
|
|
|
|
// running average
|
|
|
|
if(t>0){
|
2022-05-26 21:26:13 -04:00
|
|
|
avg_e=avg_e-(avg_e-energy)/t;
|
|
|
|
avg_a=avg_a-(avg_a-alpha)/t;
|
|
|
|
avg_en=avg_en-(avg_en-enstrophy)/t;
|
2022-05-18 09:57:06 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
if(t>0 && t%print_freq==0){
|
2022-05-26 21:26:13 -04:00
|
|
|
fprintf(stderr,"%d % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,t*delta, avg_a, avg_e, avg_en, alpha, energy, enstrophy);
|
|
|
|
printf("%8d % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,t*delta, avg_a, avg_e, avg_en, alpha, energy, enstrophy);
|
2022-05-18 09:57:06 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2022-05-27 16:09:17 -04:00
|
|
|
// save final entry to savefile
|
|
|
|
write_u(u, K1, K2, savefile);
|
|
|
|
|
2022-05-18 09:57:06 +02:00
|
|
|
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
|
|
|
|
return(0);
|
|
|
|
}
|
|
|
|
|
2022-05-18 22:22:42 +02:00
|
|
|
// compute solution as a function of time, but do not print anything (useful for debugging)
|
|
|
|
int quiet(
|
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
int N1,
|
|
|
|
int N2,
|
|
|
|
unsigned int nsteps,
|
|
|
|
double nu,
|
|
|
|
double delta,
|
2022-05-25 11:12:02 -04:00
|
|
|
double L,
|
2022-05-26 15:05:30 -04:00
|
|
|
_Complex double* u0,
|
2022-05-26 15:16:44 -04:00
|
|
|
_Complex double* g,
|
2022-05-27 16:09:17 -04:00
|
|
|
unsigned int nthreads,
|
|
|
|
FILE* savefile
|
2022-05-18 22:22:42 +02:00
|
|
|
){
|
|
|
|
_Complex double* u;
|
|
|
|
_Complex double* tmp1;
|
|
|
|
_Complex double* tmp2;
|
|
|
|
_Complex double* tmp3;
|
|
|
|
unsigned int t;
|
|
|
|
fft_vect fft1;
|
|
|
|
fft_vect fft2;
|
|
|
|
fft_vect ifft;
|
|
|
|
|
2022-05-18 23:52:01 +02:00
|
|
|
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
|
2022-05-26 15:05:30 -04:00
|
|
|
// copy initial condition
|
|
|
|
copy_u(u, u0, K1, K2);
|
2022-05-18 22:22:42 +02:00
|
|
|
|
|
|
|
// iterate
|
|
|
|
for(t=0;t<nsteps;t++){
|
2022-05-25 11:12:02 -04:00
|
|
|
ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
|
2022-05-18 22:22:42 +02:00
|
|
|
}
|
|
|
|
|
2022-05-27 16:09:17 -04:00
|
|
|
// save final entry to savefile
|
|
|
|
write_u(u, K1, K2, savefile);
|
|
|
|
|
2022-05-18 22:22:42 +02:00
|
|
|
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
|
|
|
|
return(0);
|
|
|
|
}
|
|
|
|
|
2022-05-18 09:57:06 +02:00
|
|
|
|
|
|
|
// initialize vectors for computation
|
|
|
|
int ns_init_tmps(
|
|
|
|
_Complex double ** u,
|
|
|
|
_Complex double ** tmp1,
|
|
|
|
_Complex double ** tmp2,
|
|
|
|
_Complex double ** tmp3,
|
|
|
|
fft_vect* fft1,
|
|
|
|
fft_vect* fft2,
|
|
|
|
fft_vect* ifft,
|
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
int N1,
|
2022-05-18 23:52:01 +02:00
|
|
|
int N2,
|
|
|
|
unsigned int nthreads
|
2022-05-18 09:57:06 +02:00
|
|
|
){
|
|
|
|
// velocity field
|
|
|
|
*u=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1));
|
|
|
|
|
|
|
|
// allocate tmp vectors for computation
|
|
|
|
*tmp1=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1));
|
|
|
|
*tmp2=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1));
|
|
|
|
*tmp3=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1));
|
|
|
|
|
2022-05-18 23:52:01 +02:00
|
|
|
// init threads
|
|
|
|
fftw_init_threads();
|
|
|
|
fftw_plan_with_nthreads(nthreads);
|
|
|
|
|
2022-05-18 09:57:06 +02:00
|
|
|
// prepare vectors for fft
|
|
|
|
fft1->fft=fftw_malloc(sizeof(fftw_complex)*N1*N2);
|
|
|
|
fft1->fft_plan=fftw_plan_dft_2d(N1,N2, fft1->fft, fft1->fft, FFTW_FORWARD, FFTW_MEASURE);
|
|
|
|
fft2->fft=fftw_malloc(sizeof(fftw_complex)*N1*N2);
|
|
|
|
fft2->fft_plan=fftw_plan_dft_2d(N1,N2, fft2->fft, fft2->fft, FFTW_FORWARD, FFTW_MEASURE);
|
|
|
|
ifft->fft=fftw_malloc(sizeof(fftw_complex)*N1*N2);
|
|
|
|
ifft->fft_plan=fftw_plan_dft_2d(N1,N2, ifft->fft, ifft->fft, FFTW_BACKWARD, FFTW_MEASURE);
|
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
// release vectors
|
|
|
|
int ns_free_tmps(
|
|
|
|
_Complex double* u,
|
|
|
|
_Complex double* tmp1,
|
|
|
|
_Complex double* tmp2,
|
|
|
|
_Complex double* tmp3,
|
|
|
|
fft_vect fft1,
|
|
|
|
fft_vect fft2,
|
|
|
|
fft_vect ifft
|
|
|
|
){
|
|
|
|
// free memory
|
|
|
|
fftw_destroy_plan(fft1.fft_plan);
|
|
|
|
fftw_destroy_plan(fft2.fft_plan);
|
|
|
|
fftw_destroy_plan(ifft.fft_plan);
|
|
|
|
fftw_free(fft1.fft);
|
|
|
|
fftw_free(fft2.fft);
|
|
|
|
fftw_free(ifft.fft);
|
|
|
|
|
2022-05-18 23:52:01 +02:00
|
|
|
fftw_cleanup_threads();
|
2022-05-18 09:57:06 +02:00
|
|
|
|
|
|
|
free(tmp3);
|
|
|
|
free(tmp2);
|
|
|
|
free(tmp1);
|
|
|
|
|
|
|
|
free(u);
|
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
2022-05-26 15:05:30 -04:00
|
|
|
// copy u0 to u
|
|
|
|
int copy_u(
|
2022-05-18 09:57:06 +02:00
|
|
|
_Complex double* u,
|
2022-05-26 15:05:30 -04:00
|
|
|
_Complex double* u0,
|
2022-05-18 09:57:06 +02:00
|
|
|
int K1,
|
|
|
|
int K2
|
|
|
|
){
|
2022-05-26 15:05:30 -04:00
|
|
|
int i;
|
2022-05-18 09:57:06 +02:00
|
|
|
|
2022-05-26 15:05:30 -04:00
|
|
|
for(i=0;i<(2*K1+1)*(2*K2+1);i++){
|
|
|
|
u[i]=u0[i];
|
2022-05-18 09:57:06 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2018-01-11 22:48:14 +00:00
|
|
|
// next time step for Irreversible Navier-Stokes equation
|
2022-05-18 09:57:06 +02:00
|
|
|
int ins_step(
|
|
|
|
_Complex double* u,
|
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
int N1,
|
|
|
|
int N2,
|
|
|
|
double nu,
|
|
|
|
double delta,
|
2022-05-25 11:12:02 -04:00
|
|
|
double L,
|
2022-05-26 15:16:44 -04:00
|
|
|
_Complex double* g,
|
2022-05-18 09:57:06 +02:00
|
|
|
fft_vect fft1,
|
|
|
|
fft_vect fft2,
|
|
|
|
fft_vect ifft,
|
|
|
|
_Complex double* tmp1,
|
|
|
|
_Complex double* tmp2,
|
|
|
|
_Complex double* tmp3
|
|
|
|
){
|
2018-01-11 22:48:14 +00:00
|
|
|
int kx,ky;
|
|
|
|
|
|
|
|
// k1
|
2022-05-25 11:12:02 -04:00
|
|
|
ins_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
|
2018-01-11 22:48:14 +00:00
|
|
|
// add to output
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
tmp3[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]+delta/6*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// u+h*k1/2
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
tmp2[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]+delta/2*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
// k2
|
2022-05-25 11:12:02 -04:00
|
|
|
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
|
2018-01-11 22:48:14 +00:00
|
|
|
// add to output
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
tmp3[klookup(kx,ky,2*K1+1,2*K2+1)]+=delta/3*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// u+h*k2/2
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
tmp2[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]+delta/2*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
// k3
|
2022-05-25 11:12:02 -04:00
|
|
|
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
|
2018-01-11 22:48:14 +00:00
|
|
|
// add to output
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
tmp3[klookup(kx,ky,2*K1+1,2*K2+1)]+=delta/3*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// u+h*k3
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
tmp2[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]+delta*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
// k4
|
2022-05-25 11:12:02 -04:00
|
|
|
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
|
2018-01-11 22:48:14 +00:00
|
|
|
// add to output
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
u[klookup(kx,ky,2*K1+1,2*K2+1)]=tmp3[klookup(kx,ky,2*K1+1,2*K2+1)]+delta/6*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
return(0);
|
|
|
|
}
|
|
|
|
|
|
|
|
// right side of Irreversible Navier-Stokes equation
|
2022-05-18 09:57:06 +02:00
|
|
|
int ins_rhs(
|
|
|
|
_Complex double* out,
|
|
|
|
_Complex double* u,
|
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
int N1,
|
|
|
|
int N2,
|
|
|
|
double nu,
|
2022-05-25 11:12:02 -04:00
|
|
|
double L,
|
2022-05-26 15:16:44 -04:00
|
|
|
_Complex double* g,
|
2022-05-18 09:57:06 +02:00
|
|
|
fft_vect fft1,
|
|
|
|
fft_vect fft2,
|
|
|
|
fft_vect ifft
|
|
|
|
){
|
2018-01-11 22:48:14 +00:00
|
|
|
int kx,ky;
|
2022-05-18 09:57:06 +02:00
|
|
|
int i;
|
2018-01-11 22:48:14 +00:00
|
|
|
|
2022-05-18 10:53:00 +02:00
|
|
|
// compute convolution term
|
|
|
|
ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
|
|
|
|
|
2022-05-18 22:22:42 +02:00
|
|
|
|
|
|
|
/*
|
|
|
|
// compare convolution term (store result in fft1.fft)
|
|
|
|
ns_T_nofft(fft1.fft, u, K1, K2, N1, N2);
|
|
|
|
double cmp=0.;
|
|
|
|
for(i=0;i<N1*N2;i++){
|
|
|
|
cmp+=(ifft.fft[i]-fft1.fft[i])*(ifft.fft[i]-fft1.fft[i]);
|
|
|
|
}
|
|
|
|
printf("% .15e\n",sqrt(cmp));
|
|
|
|
*/
|
|
|
|
|
|
|
|
|
2022-05-18 10:53:00 +02:00
|
|
|
for(i=0; i<(2*K1+1)*(2*K2+1); i++){
|
|
|
|
out[i]=0;
|
|
|
|
}
|
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
if(kx!=0 || ky!=0){
|
2023-03-31 15:56:50 -04:00
|
|
|
out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+g[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
|
2022-05-18 10:53:00 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
return(0);
|
|
|
|
}
|
|
|
|
|
|
|
|
// convolution term in right side of convolution equation
|
|
|
|
int ns_T(
|
|
|
|
_Complex double* u,
|
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
int N1,
|
|
|
|
int N2,
|
|
|
|
fft_vect fft1,
|
|
|
|
fft_vect fft2,
|
|
|
|
fft_vect ifft
|
|
|
|
){
|
|
|
|
int kx,ky;
|
|
|
|
int i;
|
|
|
|
|
2018-01-12 19:20:59 +00:00
|
|
|
// F(px/|p|*u)*F(qy*|q|*u)
|
2018-01-11 22:48:14 +00:00
|
|
|
// init to 0
|
2022-05-18 09:57:06 +02:00
|
|
|
for(i=0; i<N1*N2; i++){
|
|
|
|
fft1.fft[i]=0;
|
|
|
|
fft2.fft[i]=0;
|
|
|
|
ifft.fft[i]=0;
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
// fill modes
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
2018-01-11 22:48:14 +00:00
|
|
|
if(kx!=0 || ky!=0){
|
2022-05-18 10:53:00 +02:00
|
|
|
fft1.fft[klookup(kx,ky,N1,N2)]=kx/sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N1;
|
|
|
|
fft2.fft[klookup(kx,ky,N1,N2)]=ky*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2;
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2018-01-12 19:20:59 +00:00
|
|
|
|
2018-01-11 22:48:14 +00:00
|
|
|
// fft
|
2022-05-18 09:57:06 +02:00
|
|
|
fftw_execute(fft1.fft_plan);
|
|
|
|
fftw_execute(fft2.fft_plan);
|
|
|
|
// write to ifft
|
|
|
|
for(i=0;i<N1*N2;i++){
|
2022-05-25 22:32:03 -04:00
|
|
|
ifft.fft[i]=fft1.fft[i]*fft2.fft[i];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
|
2018-01-12 19:20:59 +00:00
|
|
|
// F(py/|p|*u)*F(qx*|q|*u)
|
2018-01-11 22:48:14 +00:00
|
|
|
// init to 0
|
2022-05-18 09:57:06 +02:00
|
|
|
for(i=0; i<N1*N2; i++){
|
|
|
|
fft1.fft[i]=0;
|
|
|
|
fft2.fft[i]=0;
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
// fill modes
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
2018-01-11 22:48:14 +00:00
|
|
|
if(kx!=0 || ky!=0){
|
2022-05-18 10:53:00 +02:00
|
|
|
fft1.fft[klookup(kx,ky,N1,N2)]=ky/sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N1;
|
|
|
|
fft2.fft[klookup(kx,ky,N1,N2)]=kx*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2;
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2018-01-12 19:20:59 +00:00
|
|
|
|
2018-01-11 22:48:14 +00:00
|
|
|
// fft
|
2022-05-18 09:57:06 +02:00
|
|
|
fftw_execute(fft1.fft_plan);
|
|
|
|
fftw_execute(fft2.fft_plan);
|
|
|
|
// write to ifft
|
|
|
|
for(i=0;i<N1*N2;i++){
|
2022-05-25 22:32:03 -04:00
|
|
|
ifft.fft[i]=ifft.fft[i]-fft1.fft[i]*fft2.fft[i];
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
2018-01-12 19:20:59 +00:00
|
|
|
|
2018-01-11 22:48:14 +00:00
|
|
|
// inverse fft
|
2022-05-18 09:57:06 +02:00
|
|
|
fftw_execute(ifft.fft_plan);
|
2022-05-17 14:31:22 +02:00
|
|
|
|
2023-03-31 15:56:50 -04:00
|
|
|
// enforce T(u,-k)=T(u,k)^*
|
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
ifft.fft[klookup(kx,ky,N1,N2)]=(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2018-01-11 22:48:14 +00:00
|
|
|
return(0);
|
|
|
|
}
|
|
|
|
|
2022-05-18 22:22:42 +02:00
|
|
|
// convolution term in right side of convolution equation, computed without fourier transform
|
|
|
|
int ns_T_nofft(
|
|
|
|
_Complex double* out,
|
|
|
|
_Complex double* u,
|
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
int N1,
|
|
|
|
int N2
|
|
|
|
){
|
|
|
|
int kx,ky;
|
|
|
|
int px,py;
|
|
|
|
int qx,qy;
|
|
|
|
|
|
|
|
// loop over K's (needs N1>=4*K1+1 and N2>=4*K2+1)
|
|
|
|
if (N1<4*K1+1 || N2<4*K2+1){
|
|
|
|
fprintf(stderr,"error: N1 and N2 need t be >= 4*K1+1 and 4*K2+1 respectively\n");
|
|
|
|
return(-1);
|
|
|
|
}
|
|
|
|
for(kx=-2*K1;kx<=2*K1;kx++){
|
|
|
|
for(ky=-2*K2;ky<=2*K2;ky++){
|
|
|
|
// init
|
|
|
|
out[klookup(kx,ky,N1,N2)]=0.;
|
|
|
|
|
|
|
|
for(px=-K1;px<=K1;px++){
|
|
|
|
for(py=-K2;py<=K2;py++){
|
|
|
|
qx=kx-px;
|
|
|
|
qy=ky-py;
|
|
|
|
|
|
|
|
// cutoff in q
|
|
|
|
if(qx>=-K1 && qx<=K1 && qy>=-K2 && qy<=K2 && qx*qx+qy*qy>0 && px*px+py*py>0){
|
|
|
|
out[klookup(kx,ky,N1,N2)]+=(-qx*py+qy*px)*sqrt(qx*qx+qy*qy)/sqrt(px*px+py*py)*u[klookup(px,py,2*K1+1,2*K2+1)]*u[klookup(qx,qy,2*K1+1,2*K2+1)];
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2023-03-31 15:56:50 -04:00
|
|
|
// enforce T(u,-k)=T(u,k)^*
|
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
out[klookup(kx,ky,N1,N2)]=(out[klookup(kx,ky,N1,N2)]+conj(out[klookup(-kx,-ky,N1,N2)]))/2;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2022-05-18 22:22:42 +02:00
|
|
|
return 0;
|
|
|
|
}
|
2018-01-11 22:48:14 +00:00
|
|
|
|
|
|
|
// compute alpha
|
2022-05-26 21:26:13 -04:00
|
|
|
double compute_alpha(
|
2022-05-18 09:57:06 +02:00
|
|
|
_Complex double* u,
|
|
|
|
int K1,
|
|
|
|
int K2,
|
2022-05-26 15:16:44 -04:00
|
|
|
_Complex double* g
|
2022-05-18 09:57:06 +02:00
|
|
|
){
|
2018-01-11 22:48:14 +00:00
|
|
|
_Complex double num=0;
|
|
|
|
_Complex double denom=0;
|
|
|
|
int kx,ky;
|
|
|
|
|
2022-05-18 09:57:06 +02:00
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for(ky=-K2;ky<=K2;ky++){
|
|
|
|
denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)])*(1+(ky!=0?kx*kx/ky/ky:0));
|
2022-05-26 15:16:44 -04:00
|
|
|
num+=(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(g[klookup(kx,ky,2*K1+1,2*K2+1)])*(1+(ky!=0?kx*kx/ky/ky:0));
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2022-05-26 21:26:13 -04:00
|
|
|
return __real__ (num/denom);
|
|
|
|
}
|
|
|
|
|
|
|
|
// compute energy
|
|
|
|
double compute_energy(
|
|
|
|
_Complex double* u,
|
|
|
|
int K1,
|
|
|
|
int K2
|
|
|
|
){
|
|
|
|
int kx,ky;
|
|
|
|
double out=0.;
|
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for (ky=-K2;ky<=K2;ky++){
|
|
|
|
out+=__real__ (u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]));
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return out;
|
|
|
|
}
|
|
|
|
|
|
|
|
// compute enstrophy
|
|
|
|
double compute_enstrophy(
|
|
|
|
_Complex double* u,
|
|
|
|
int K1,
|
|
|
|
int K2,
|
|
|
|
double L
|
|
|
|
){
|
|
|
|
int kx,ky;
|
|
|
|
double out=0.;
|
|
|
|
for(kx=-K1;kx<=K1;kx++){
|
|
|
|
for (ky=-K2;ky<=K2;ky++){
|
|
|
|
out+=__real__ (4*M_PI*M_PI/L/L*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]));
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return out;
|
2018-01-11 22:48:14 +00:00
|
|
|
}
|
2022-05-18 09:57:06 +02:00
|
|
|
|
|
|
|
|
|
|
|
// get index for kx,ky in array of size S
|
|
|
|
int klookup(
|
|
|
|
int kx,
|
|
|
|
int ky,
|
|
|
|
int S1,
|
|
|
|
int S2
|
|
|
|
){
|
|
|
|
return (kx>=0 ? kx : S1+kx)*S2 + (ky>=0 ? ky : S2+ky);
|
|
|
|
}
|