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

add L parameter

Browse Source
This commit is contained in:
Ian Jauslin 2022-05-25 11:12:02 -04:00
parent d37d6104e0
commit a55065f474
3 changed files with 43 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
src/main.c
View File

@ -14,8 +14,8 @@ int print_usage();
// read command line arguments // read command line arguments
int read_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* nthreads); int read_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* nthreads);
int read_params(char* params, int* K1, int* K2, int* N1, int* N2, unsigned int* nsteps, double* nu, double* delta, unsigned int* print_freq); int read_params(char* params, int* K1, int* K2, int* N1, int* N2, unsigned int* nsteps, double* nu, double* delta, double* L, unsigned int* print_freq);
int set_parameter(char* lhs, char* rhs, int* K1, int* K2, int* N1, int* N2, unsigned int* nsteps, double* nu, double* delta, unsigned int* print_freq, bool* setN1, bool* setN2); int set_parameter(char* lhs, char* rhs, int* K1, int* K2, int* N1, int* N2, unsigned int* nsteps, double* nu, double* delta, double* L, unsigned int* print_freq, bool* setN1, bool* setN2);
#define COMMAND_UK 1 #define COMMAND_UK 1
@ -35,7 +35,7 @@ int main (
int K1,K2; int K1,K2;
int N1,N2; int N1,N2;
unsigned int nsteps; unsigned int nsteps;
double nu,delta; double nu,delta,L;
_Complex double (*g)(int,int); _Complex double (*g)(int,int);
int ret; int ret;
unsigned int driving,command; unsigned int driving,command;
@ -51,7 +51,7 @@ int main (
return(-1); return(-1);
} }
// read params // read params
ret=read_params(params, &K1, &K2, &N1, &N2, &nsteps, &nu, &delta, &print_freq); ret=read_params(params, &K1, &K2, &N1, &N2, &nsteps, &nu, &delta, &L, &print_freq);
if(ret<0){ if(ret<0){
return(-1); return(-1);
} }
@ -71,16 +71,16 @@ int main (
// run command // run command
if (command==COMMAND_UK){ if (command==COMMAND_UK){
uk(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq, nthreads); uk(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
} }
else if (command==COMMAND_ENERGY){ else if (command==COMMAND_ENERGY){
energy(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq, nthreads); energy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
} }
else if(command==COMMAND_ENSTROPHY){ else if(command==COMMAND_ENSTROPHY){
energy(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq, nthreads); energy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
} }
else if(command==COMMAND_QUIET){ else if(command==COMMAND_QUIET){
quiet(K1, K2, N1, N2, nsteps, nu, delta, g, nthreads); quiet(K1, K2, N1, N2, nsteps, nu, delta, L, g, nthreads);
} }
else if(command==0){ else if(command==0){
fprintf(stderr, "error: no command specified\n"); fprintf(stderr, "error: no command specified\n");
@ -201,6 +201,7 @@ int read_params(
unsigned int* nsteps, unsigned int* nsteps,
double* nu, double* nu,
double* delta, double* delta,
double* L,
unsigned int* print_freq unsigned int* print_freq
){ ){
int ret; int ret;
@ -222,6 +223,7 @@ int read_params(
*delta=0.0001220703125; *delta=0.0001220703125;
//nu=2^-11 //nu=2^-11
*nu=0.00048828125; *nu=0.00048828125;
*L=2*M_PI;
*nsteps=10000000; *nsteps=10000000;
*print_freq=1000; *print_freq=1000;
@ -244,7 +246,7 @@ int read_params(
break; break;
case ';': case ';':
//set parameter //set parameter
ret=set_parameter(buffer_lhs,buffer_rhs,K1,K2,N1,N2,nsteps,nu,delta,print_freq,&setN1,&setN2); ret=set_parameter(buffer_lhs,buffer_rhs,K1,K2,N1,N2,nsteps,nu,delta,L,print_freq,&setN1,&setN2);
if(ret<0){ if(ret<0){
return ret; return ret;
} }
@ -271,7 +273,7 @@ int read_params(
// set last param // set last param
if (*params!='\0'){ if (*params!='\0'){
ret=set_parameter(buffer_lhs,buffer_rhs,K1,K2,N1,N2,nsteps,nu,delta,print_freq,&setN1,&setN2); ret=set_parameter(buffer_lhs,buffer_rhs,K1,K2,N1,N2,nsteps,nu,delta,L,print_freq,&setN1,&setN2);
if(ret<0){ if(ret<0){
return ret; return ret;
} }
@ -305,6 +307,7 @@ int set_parameter(
unsigned int* nsteps, unsigned int* nsteps,
double* nu, double* nu,
double* delta, double* delta,
double* L,
unsigned int* print_freq, unsigned int* print_freq,
bool* setN1, bool* setN1,
bool* setN2 bool* setN2
@ -380,6 +383,13 @@ int set_parameter(
return(-1); return(-1);
} }
} }
else if (strcmp(lhs,"L")==0){
ret=sscanf(rhs,"%lf",L);
if(ret!=1){
fprintf(stderr, "error: parameter 'L' should be a double\n got '%s'\n",rhs);
return(-1);
}
}
else if (strcmp(lhs,"print_freq")==0){ else if (strcmp(lhs,"print_freq")==0){
ret=sscanf(rhs,"%u",print_freq); ret=sscanf(rhs,"%u",print_freq);
if(ret!=1){ if(ret!=1){

26
src/navier-stokes.c
View File

@ -2,8 +2,6 @@
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
#define M_PI 3.14159265358979323846
// compute solution as a function of time // compute solution as a function of time
int uk( int uk(
int K1, int K1,
@ -13,6 +11,7 @@ int uk(
unsigned int nsteps, unsigned int nsteps,
double nu, double nu,
double delta, double delta,
double L,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
unsigned int print_freq, unsigned int print_freq,
unsigned int nthreads unsigned int nthreads
@ -44,7 +43,7 @@ int uk(
// iterate // iterate
for(t=0;t<nsteps;t++){ for(t=0;t<nsteps;t++){
ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
if(t%print_freq==0){ if(t%print_freq==0){
fprintf(stderr,"%d % .8e ",t,t*delta); fprintf(stderr,"%d % .8e ",t,t*delta);
@ -77,6 +76,7 @@ int energy(
unsigned int nsteps, unsigned int nsteps,
double nu, double nu,
double delta, double delta,
double L,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
unsigned int print_freq, unsigned int print_freq,
unsigned int nthreads unsigned int nthreads
@ -97,7 +97,7 @@ int energy(
// iterate // iterate
for(t=0;t<nsteps;t++){ for(t=0;t<nsteps;t++){
ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
if(t%print_freq==0){ if(t%print_freq==0){
@ -126,6 +126,7 @@ int enstrophy(
unsigned int nsteps, unsigned int nsteps,
double nu, double nu,
double delta, double delta,
double L,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
unsigned int print_freq, unsigned int print_freq,
unsigned int nthreads unsigned int nthreads
@ -150,7 +151,7 @@ int enstrophy(
// iterate // iterate
for(t=0;t<nsteps;t++){ for(t=0;t<nsteps;t++){
ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
alpha=compute_alpha(u, K1, K2, g); alpha=compute_alpha(u, K1, K2, g);
// running average // running average
@ -177,6 +178,7 @@ int quiet(
unsigned int nsteps, unsigned int nsteps,
double nu, double nu,
double delta, double delta,
double L,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
unsigned int nthreads unsigned int nthreads
){ ){
@ -194,7 +196,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, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
} }
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft); ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
@ -342,6 +344,7 @@ int ins_step(
int N2, int N2,
double nu, double nu,
double delta, double delta,
double L,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
fft_vect fft1, fft_vect fft1,
fft_vect fft2, fft_vect fft2,
@ -353,7 +356,7 @@ int ins_step(
int kx,ky; int kx,ky;
// k1 // k1
ins_rhs(tmp1, u, K1, K2, N1, N2, nu, g, fft1, fft2, ifft); ins_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
// 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++){
@ -368,7 +371,7 @@ int ins_step(
} }
} }
// k2 // k2
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft); ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
// 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++){
@ -383,7 +386,7 @@ int ins_step(
} }
} }
// k3 // k3
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft); ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
// 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++){
@ -398,7 +401,7 @@ int ins_step(
} }
} }
// k4 // k4
ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft); ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
// 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++){
@ -418,6 +421,7 @@ int ins_rhs(
int N1, int N1,
int N2, int N2,
double nu, double nu,
double L,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
fft_vect fft1, fft_vect fft1,
fft_vect fft2, fft_vect fft2,
@ -447,7 +451,7 @@ int ins_rhs(
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*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/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*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
} }
} }
} }

15
src/navier-stokes.h
View File

@ -4,6 +4,7 @@
#include <complex.h> #include <complex.h>
#include <fftw3.h> #include <fftw3.h>
#define M_PI 3.14159265358979323846
// arrays on which the ffts are performed // arrays on which the ffts are performed
typedef struct fft_vects { typedef struct fft_vects {
@ -12,16 +13,16 @@ 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, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads); int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
// compute the energy as a function of time // compute the energy as a function of time
int energy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads); int energy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
// compute enstrophy // compute enstrophy
int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads); int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
// 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, _Complex double (*g)(int,int), unsigned int nthreads); int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double (*g)(int,int), unsigned int nthreads);
// initialize vectors for computation // initialize vectors for computation
@ -33,15 +34,15 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
int ns_init_u( _Complex double* u, int K1, int K2); int ns_init_u( _Complex double* u, int K1, int K2);
// next time step for Irreversible Navier-Stokes equation // next time step for Irreversible Navier-Stokes equation
int ins_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3); int ins_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3);
// right side of Irreversible Navier-Stokes equation // right side of Irreversible Navier-Stokes equation
int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2, fft_vect ifft); int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2, fft_vect ifft);
// 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);
// convolution term in right side of equation // convolution term in right side of equation (computed without 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