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

choose initial condition on cli

Browse Source
This commit is contained in:
Ian Jauslin 2022-05-26 15:05:30 -04:00
parent 8877b63549
commit d4254c6b8e
5 changed files with 172 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
src/init.c Normal file
View File

@ -0,0 +1,62 @@
#include "init.h"
#include "navier-stokes.h"
#include <stdlib.h>
#include <math.h>
// random initial condition
int init_random (
_Complex double* u0,
int K1,
int K2,
int seed
){
int kx,ky;
double rescale;
double x,y;
srand(seed);
// random init (set half, then the other half are the conjugates)
for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
if (kx!=0 || ky>0){
x=-0.5+((double) rand())/RAND_MAX;
y=-0.5+((double) rand())/RAND_MAX;
u0[klookup(kx,ky,2*K1+1,2*K2+1)]=x+y*I;
u0[klookup(-kx,-ky,2*K1+1,2*K2+1)]=conj(u0[klookup(kx,ky,2*K1+1,2*K2+1)]);
}
}
}
// rescale to match with Gallavotti's initialization
rescale=0;
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
rescale=rescale+((__real__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])*(__real__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])+(__imag__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])*(__imag__ u0[klookup(kx,ky,2*K1+1,2*K2+1)]))*(kx*kx+ky*ky);
}
}
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
u0[klookup(kx,ky,2*K1+1,2*K2+1)]=u0[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(1.54511597324389e+02/rescale);
}
}
return 0;
}
// Gaussian initial condition
int init_gaussian (
_Complex double* u0,
int K1,
int K2
){
int kx,ky;
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
u0[klookup(kx,ky,2*K1+1,2*K2+1)]=(kx*kx+ky*ky)*exp(-(kx*kx+ky*ky));
}
}
return 0;
}

10
src/init.h Normal file
View File

@ -0,0 +1,10 @@
#ifndef INIT_H
#define INIT_H
// random initial condition
int init_random(_Complex double* u0, int K1, int K2, int seed);
// Gaussian initial condition
int init_gaussian(_Complex double* u0, int K1, int K2);
#endif

84
src/main.c
View File

@ -8,6 +8,7 @@
#include <stdbool.h> #include <stdbool.h>
#include "navier-stokes.h" #include "navier-stokes.h"
#include "driving.h" #include "driving.h"
#include "init.h"
// 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 {
@ -20,16 +21,20 @@ typedef struct nstrophy_parameters {
double delta; double delta;
double L; double L;
unsigned int print_freq; unsigned int print_freq;
int seed;
} nstrophy_parameters; } nstrophy_parameters;
// usage message // usage message
int print_usage(); 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* init, unsigned int* nthreads);
int read_params(char* param_str, nstrophy_parameters* parameters); int read_params(char* param_str, nstrophy_parameters* parameters);
int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2); int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2);
// set initial condition
_Complex double* set_init(unsigned int init, nstrophy_parameters parameters);
#define COMMAND_UK 1 #define COMMAND_UK 1
#define COMMAND_ENSTROPHY 2 #define COMMAND_ENSTROPHY 2
#define COMMAND_QUIET 3 #define COMMAND_QUIET 3
@ -38,6 +43,9 @@ int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* s
#define DRIVING_ZERO 1 #define DRIVING_ZERO 1
#define DRIVING_TEST 2 #define DRIVING_TEST 2
#define INIT_RANDOM 1
#define INIT_GAUSSIAN 2
int main ( int main (
int argc, int argc,
@ -47,22 +55,26 @@ int main (
nstrophy_parameters parameters; nstrophy_parameters parameters;
_Complex double (*g)(int,int); _Complex double (*g)(int,int);
int ret; int ret;
unsigned int driving,command; unsigned int driving,command,init;
unsigned int nthreads=1; unsigned int nthreads=1;
_Complex double* u0;
command=0; command=0;
driving=0; driving=0;
init=0;
// read command line arguments // read command line arguments
ret=read_args(argc, argv, &param_str, &driving, &command, &nthreads); ret=read_args(argc, argv, &param_str, &driving, &init, &command, &nthreads);
if(ret<0){ if(ret<0){
return(-1); return(-1);
} }
// read params // read params
ret=read_params(param_str, &parameters); ret=read_params(param_str, &parameters);
if(ret<0){ if(ret<0){
return(-1); return(-1);
} }
// set driving force // set driving force
switch(driving){ switch(driving){
case DRIVING_ZERO: case DRIVING_ZERO:
@ -77,30 +89,35 @@ int main (
break; break;
} }
// set initial condition
u0=set_init(init, parameters);
// 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, g, parameters.print_freq, nthreads); uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads);
} }
else if (command==COMMAND_ENERGY){ else if (command==COMMAND_ENERGY){
energy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, g, parameters.print_freq, nthreads); energy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads);
} }
else if(command==COMMAND_ENSTROPHY){ else if(command==COMMAND_ENSTROPHY){
enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, g, parameters.print_freq, nthreads); enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads);
} }
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, g, nthreads); quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, nthreads);
} }
else if(command==0){ else if(command==0){
fprintf(stderr, "error: no command specified\n"); fprintf(stderr, "error: no command specified\n");
print_usage(); print_usage();
} }
free(u0);
return(0); return(0);
} }
// usage message // usage message
int print_usage(){ int print_usage(){
fprintf(stderr, "usage:\n nstrophy [-p parameters] [-g driving_force] <command>\n\n"); fprintf(stderr, "usage:\n nstrophy [-t nthreads] [-p parameters] [-g driving_force] [-i initial_condition] <command>\n\n");
return(0); return(0);
} }
@ -108,11 +125,13 @@ int print_usage(){
#define CP_FLAG_PARAMS 1 #define CP_FLAG_PARAMS 1
#define CP_FLAG_DRIVING 2 #define CP_FLAG_DRIVING 2
#define CP_FLAG_NTHREADS 3 #define CP_FLAG_NTHREADS 3
#define CP_FLAG_INIT 4
int read_args( int read_args(
int argc, int argc,
const char* argv[], const char* argv[],
char** params, char** params,
unsigned int* driving_force, unsigned int* driving_force,
unsigned int* init,
unsigned int* command, unsigned int* command,
unsigned int* nthreads unsigned int* nthreads
){ ){
@ -138,6 +157,9 @@ int read_args(
case 't': case 't':
flag=CP_FLAG_NTHREADS; flag=CP_FLAG_NTHREADS;
break; break;
case 'i':
flag=CP_FLAG_INIT;
break;
default: default:
fprintf(stderr, "unrecognized option '-%c'\n", *ptr); fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
print_usage(); print_usage();
@ -174,6 +196,20 @@ int read_args(
} }
flag=0; flag=0;
} }
// initial condition
else if(flag==CP_FLAG_INIT){
if (strcmp(argv[i],"random")==0){
*init=INIT_RANDOM;
}
else if (strcmp(argv[i],"gaussian")==0){
*init=INIT_GAUSSIAN;
}
else{
fprintf(stderr, "error: unrecognized initial condition '%s'\n",argv[i]);
return(-1);
}
flag=0;
}
// computation to run // computation to run
else{ else{
if(strcmp(argv[i], "uk")==0){ if(strcmp(argv[i], "uk")==0){
@ -226,6 +262,7 @@ int read_params(
parameters->L=2*M_PI; parameters->L=2*M_PI;
parameters->nsteps=10000000; parameters->nsteps=10000000;
parameters->print_freq=1000; parameters->print_freq=1000;
parameters->seed=17;
if (param_str!=NULL){ if (param_str!=NULL){
// init // init
@ -389,6 +426,13 @@ int set_parameter(
return(-1); return(-1);
} }
} }
else if (strcmp(lhs,"random_seed")==0){
ret=sscanf(rhs,"%d",&(parameters->seed));
if(ret!=1){
fprintf(stderr, "error: parameter 'random_seed' should be an integer\n got '%s'\n",rhs);
return(-1);
}
}
else{ else{
fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs); fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs);
return(-1); return(-1);
@ -396,3 +440,27 @@ int set_parameter(
return(0); return(0);
} }
// set initial condition
_Complex double* set_init(
unsigned int init,
nstrophy_parameters parameters
){
_Complex double* u0=calloc(sizeof(_Complex double),(2*parameters.K1+1)*(2*parameters.K2+1));
switch(init){
case INIT_RANDOM:
init_random(u0, parameters.K1, parameters.K2, parameters.seed);
break;
case INIT_GAUSSIAN:
init_gaussian(u0, parameters.K1, parameters.K2);
break;
default:
init_gaussian(u0, parameters.K1, parameters.K2);
break;
}
return u0;
}

70
src/navier-stokes.c
View File

@ -12,6 +12,7 @@ int uk(
double nu, double nu,
double delta, double delta,
double L, double L,
_Complex double* u0,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
unsigned int print_freq, unsigned int print_freq,
unsigned int nthreads unsigned int nthreads
@ -27,7 +28,8 @@ int uk(
int kx,ky; int kx,ky;
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads); ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
ns_init_u(u, K1, K2); // copy initial condition
copy_u(u, u0, K1, K2);
// print column headers // print column headers
printf("# 1:i 2:t "); printf("# 1:i 2:t ");
@ -77,6 +79,7 @@ int energy(
double nu, double nu,
double delta, double delta,
double L, double L,
_Complex double* u0,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
unsigned int print_freq, unsigned int print_freq,
unsigned int nthreads unsigned int nthreads
@ -93,7 +96,8 @@ int energy(
double energy; double energy;
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads); ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
ns_init_u(u, K1, K2); // copy initial condition
copy_u(u, u0, K1, K2);
// iterate // iterate
for(t=0;t<nsteps;t++){ for(t=0;t<nsteps;t++){
@ -126,6 +130,7 @@ int enstrophy(
double nu, double nu,
double delta, double delta,
double L, double L,
_Complex double* u0,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
unsigned int print_freq, unsigned int print_freq,
unsigned int nthreads unsigned int nthreads
@ -142,7 +147,8 @@ int enstrophy(
fft_vect ifft; fft_vect ifft;
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads); ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
ns_init_u(u, K1, K2); // copy initial condition
copy_u(u, u0, K1, K2);
// init running average // init running average
@ -178,6 +184,7 @@ int quiet(
double nu, double nu,
double delta, double delta,
double L, double L,
_Complex double* u0,
_Complex double (*g)(int,int), _Complex double (*g)(int,int),
unsigned int nthreads unsigned int nthreads
){ ){
@ -191,7 +198,8 @@ int quiet(
fft_vect ifft; fft_vect ifft;
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads); ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
ns_init_u(u, K1, K2); // copy initial condition
copy_u(u, u0, K1, K2);
// iterate // iterate
for(t=0;t<nsteps;t++){ for(t=0;t<nsteps;t++){
@ -272,59 +280,17 @@ int ns_free_tmps(
// initial value // copy u0 to u
int ns_init_u( int copy_u(
_Complex double* u, _Complex double* u,
_Complex double* u0,
int K1, int K1,
int K2 int K2
){ ){
int kx,ky; int i;
/* for(i=0;i<(2*K1+1)*(2*K2+1);i++){
srand(17); u[i]=u0[i];
// random init (set half, then the other half are the conjugates)
for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
if (kx!=0 || ky>0){
double x=-0.5+((double) rand())/RAND_MAX;
double y=-0.5+((double) rand())/RAND_MAX;
u[klookup(kx,ky,2*K1+1,2*K2+1)]=x+y*I;
u[klookup(-kx,-ky,2*K1+1,2*K2+1)]=conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]);
}
}
}
// rescale to match with Gallavotti's initialization
double rescale;
rescale=0;
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
rescale=rescale+((__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)])*(__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)])+(__imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)])*(__imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]))*(kx*kx+ky*ky);
}
}
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
u[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(1.54511597324389e+02/rescale);
}
}
*/
/*
// constant init
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
u[klookup(kx,ky,2*K1+1,2*K2+1)]=1.;
}
}
*/
// gaussian init
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
u[klookup(kx,ky,2*K1+1,2*K2+1)]=(kx*kx+ky*ky)*exp(-(kx*kx+ky*ky));
}
} }
return 0; return 0;

12
src/navier-stokes.h
View File

@ -13,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, double L, _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* u0, _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, double L, _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* u0, _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, double L, _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* u0, _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, double L, _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* u0, _Complex double (*g)(int,int), unsigned int nthreads);
// initialize vectors for computation // initialize vectors for computation
@ -30,8 +30,8 @@ int ns_init_tmps( _Complex double **u, _Complex double ** tmp1, _Complex double
// release vectors // 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); 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);
// initial value // copy u0 to u
int ns_init_u( _Complex double* u, 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 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)(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);