software/Nstrophy
1
Fork 0
Code Issues Pull Requests Releases Activity

Reversible equation

Browse Source
This commit is contained in:
Ian Jauslin 2023-04-05 20:33:38 -04:00
parent 75de7e03b7
commit 0bf223bcb9
3 changed files with 87 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/main.c
View File

@ -14,6 +14,7 @@
// structure to store parameters, to make it easier to pass parameters to CLI functions // structure to store parameters, to make it easier to pass parameters to CLI functions
typedef struct nstrophy_parameters { typedef struct nstrophy_parameters {
bool irreversible;
int K1; int K1;
int K2; int K2;
int N1; int N1;
@ -120,13 +121,13 @@ int main (
// run command // run command
if (command==COMMAND_UK){ if (command==COMMAND_UK){
uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, parameters.starting_time, nthreads, savefile); uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.starting_time, nthreads, savefile);
} }
else if(command==COMMAND_EEA){ else if(command==COMMAND_EEA){
eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, parameters.starting_time, nthreads, savefile); eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.starting_time, nthreads, savefile);
} }
else if(command==COMMAND_QUIET){ else if(command==COMMAND_QUIET){
quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, nthreads, savefile); quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, nthreads, savefile);
} }
else if(command==0){ else if(command==0){
fprintf(stderr, "error: no command specified\n"); fprintf(stderr, "error: no command specified\n");
@ -337,6 +338,7 @@ int read_params(
bool lhs=true; bool lhs=true;
// defaults // defaults
parameters->irreversible=true;
parameters->K1=16; parameters->K1=16;
parameters->K2=parameters->K1; parameters->K2=parameters->K1;
//delta=2^-13 //delta=2^-13
@ -428,7 +430,19 @@ int set_parameter(
){ ){
int ret; int ret;
if (strcmp(lhs,"K1")==0){ if (strcmp(lhs,"equation")==0){
if (strcmp(rhs,"irreversible")==0){
parameters->irreversible=true;
}
else if (strcmp(rhs,"reversible")==0){
parameters->irreversible=false;
}
else {
fprintf(stderr, "error: 'equation' should be 'irreversible' or 'reversible'\n got '%s'\n",rhs);
return(-1);
}
}
else if (strcmp(lhs,"K1")==0){
ret=sscanf(rhs,"%d",&(parameters->K1)); ret=sscanf(rhs,"%d",&(parameters->K1));
if(ret!=1){ if(ret!=1){
fprintf(stderr, "error: parameter 'K1' should be an integer\n got '%s'\n",rhs); fprintf(stderr, "error: parameter 'K1' should be an integer\n got '%s'\n",rhs);

77
src/navier-stokes.c
View File

@ -15,6 +15,7 @@ int uk(
double L, double L,
_Complex double* u0, _Complex double* u0,
_Complex double* g, _Complex double* g,
bool irreversible,
unsigned int print_freq, unsigned int print_freq,
unsigned int starting_time, unsigned int starting_time,
unsigned int nthreads, unsigned int nthreads,
@ -48,7 +49,7 @@ int uk(
// iterate // iterate
for(t=starting_time;t<starting_time+nsteps;t++){ for(t=starting_time;t<starting_time+nsteps;t++){
ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
if(t%print_freq==0){ if(t%print_freq==0){
fprintf(stderr,"%d % .8e ",t,t*delta); fprintf(stderr,"%d % .8e ",t,t*delta);
@ -87,6 +88,7 @@ int eea(
double L, double L,
_Complex double* u0, _Complex double* u0,
_Complex double* g, _Complex double* g,
bool irreversible,
unsigned int print_freq, unsigned int print_freq,
unsigned int starting_time, unsigned int starting_time,
unsigned int nthreads, unsigned int nthreads,
@ -115,10 +117,13 @@ int eea(
// iterate // iterate
for(t=starting_time;t<starting_time+nsteps;t++){ for(t=starting_time;t<starting_time+nsteps;t++){
ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
// convolution term
ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
energy=compute_energy(u, K1, K2); energy=compute_energy(u, K1, K2);
alpha=compute_alpha(u, K1, K2, g); alpha=compute_alpha(u, K1, K2, N1, N2, g, L, ifft.fft);
enstrophy=compute_enstrophy(u, K1, K2, L); enstrophy=compute_enstrophy(u, K1, K2, L);
// running average // running average
@ -153,6 +158,7 @@ int quiet(
double L, double L,
_Complex double* u0, _Complex double* u0,
_Complex double* g, _Complex double* g,
bool irreversible,
unsigned int nthreads, unsigned int nthreads,
FILE* savefile FILE* savefile
){ ){
@ -171,7 +177,7 @@ int quiet(
// iterate // iterate
for(t=0;t<nsteps;t++){ for(t=0;t<nsteps;t++){
ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
} }
// save final entry to savefile // save final entry to savefile
@ -267,8 +273,8 @@ int copy_u(
return 0; return 0;
} }
// next time step for Irreversible Navier-Stokes equation // next time step
int ins_step( int ns_step(
_Complex double* u, _Complex double* u,
int K1, int K1,
int K2, int K2,
@ -283,12 +289,13 @@ int ins_step(
fft_vect ifft, fft_vect ifft,
_Complex double* tmp1, _Complex double* tmp1,
_Complex double* tmp2, _Complex double* tmp2,
_Complex double* tmp3 _Complex double* tmp3,
bool irreversible
){ ){
int kx,ky; int kx,ky;
// k1 // k1
ins_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft); ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
@ -303,7 +310,7 @@ int ins_step(
} }
} }
// k2 // k2
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
@ -318,7 +325,7 @@ int ins_step(
} }
} }
// k3 // k3
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
@ -333,7 +340,7 @@ int ins_step(
} }
} }
// k4 // k4
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
@ -344,8 +351,8 @@ int ins_step(
return(0); return(0);
} }
// right side of Irreversible Navier-Stokes equation // right side of Irreversible/Reversible Navier-Stokes equation
int ins_rhs( int ns_rhs(
_Complex double* out, _Complex double* out,
_Complex double* u, _Complex double* u,
int K1, int K1,
@ -357,14 +364,21 @@ int ins_rhs(
_Complex double* g, _Complex double* g,
fft_vect fft1, fft_vect fft1,
fft_vect fft2, fft_vect fft2,
fft_vect ifft fft_vect ifft,
bool irreversible
){ ){
int kx,ky; int kx,ky;
int i; int i;
double alpha;
// compute convolution term // compute convolution term
ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft); ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
if (irreversible) {
alpha=nu;
} else {
alpha=compute_alpha(u,K1,K2,N1,N2,g,L,ifft.fft);
}
/* /*
// compare convolution term (store result in fft1.fft) // compare convolution term (store result in fft1.fft)
@ -376,14 +390,13 @@ int ins_rhs(
printf("% .15e\n",sqrt(cmp)); printf("% .15e\n",sqrt(cmp));
*/ */
for(i=0; i<(2*K1+1)*(2*K2+1); i++){ for(i=0; i<(2*K1+1)*(2*K2+1); i++){
out[i]=0; out[i]=0;
} }
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
if(kx!=0 || ky!=0){ if(kx!=0 || ky!=0){
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)]; out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*alpha*(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)];
} }
} }
} }
@ -514,6 +527,36 @@ int ns_T_nofft(
return 0; return 0;
} }
// compute alpha
double compute_alpha(
_Complex double* u,
int K1,
int K2,
int N1,
int N2,
_Complex double* g,
double L,
_Complex double* T
){
_Complex double num=0;
double denom=0;
int kx,ky;
num=0.;
denom=0.;
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
num+=(L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*g[klookup(kx,ky,2*K1+1,2*K2+1)]+sqrt(kx*kx+ky*ky)*T[klookup(kx,ky,N1,N2)])*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]);
denom+=__real__ (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)]);
}
}
return __real__ num/denom;
}
/*
// compute alpha // compute alpha
double compute_alpha( double compute_alpha(
_Complex double* u, _Complex double* u,
@ -533,7 +576,7 @@ double compute_alpha(
} }
return __real__ (num/denom); return __real__ (num/denom);
} }*/
// compute energy // compute energy
double compute_energy( double compute_energy(

17
src/navier-stokes.h
View File

@ -2,6 +2,7 @@
#define NAVIERSTOKES_H #define NAVIERSTOKES_H
#include <complex.h> #include <complex.h>
#include <stdbool.h>
#include <fftw3.h> #include <fftw3.h>
#define M_PI 3.14159265358979323846 #define M_PI 3.14159265358979323846
@ -13,13 +14,13 @@ typedef struct fft_vects {
} fft_vect; } fft_vect;
// compute u_k // compute u_k
int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int starting_time, unsigned int nthreadsl, FILE* savefile); int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int print_freq, unsigned int starting_time, unsigned int nthreadsl, FILE* savefile);
// compute energy, enstrophy and alpha // compute energy, enstrophy and alpha
int eea( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int starting_time, unsigned int nthreads, FILE* savefile); int eea( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int print_freq, unsigned int starting_time, unsigned int nthreads, FILE* savefile);
// compute solution as a function of time, but do not print anything (useful for debugging) // 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, double L, _Complex double* u0, _Complex double* g, unsigned int nthreads, FILE* savefile); int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int nthreads, FILE* savefile);
// initialize vectors for computation // initialize vectors for computation
@ -30,11 +31,11 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
// copy u0 to u // copy u0 to u
int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2); int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);
// next time step for Irreversible Navier-Stokes equation // next time step for Irreversible/reversible Navier-Stokes equation
int ins_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3); int ns_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3, bool irreversible);
// right side of Irreversible Navier-Stokes equation // right side of Irreversible/reversible Navier-Stokes equation
int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft); int ns_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, bool irreversible);
// convolution term in right side of equation // convolution term in right side of 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 ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft_vect fft2, fft_vect ifft);
@ -43,7 +44,7 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2); int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
// compute alpha // compute alpha
double compute_alpha( _Complex double* u, int K1, int K2, _Complex double* g); double compute_alpha( _Complex double* u, int K1, int K2, int N1, int N2, _Complex double* g, double L, _Complex double* T);
// compute energy // compute energy
double compute_energy( _Complex double* u, int K1, int K2); double compute_energy( _Complex double* u, int K1, int K2);
// compute enstrophy // compute enstrophy