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

savefile and initfile

Browse Source
This commit is contained in:
Ian Jauslin
2022-05-27 16:09:17 -04:00
parent 6b96c2d946
commit 0cdb914b57
7 changed files with 231 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/init.c
View File

@@ -1,5 +1,6 @@
#include "init.h"
#include "navier-stokes.h"
#include "io.h"
#include <stdlib.h>
#include <math.h>
@@ -60,3 +61,14 @@ int init_gaussian (
return 0;
}
// Initialize from file
int init_file (
_Complex double* u0,
int K1,
int K2,
FILE* initfile
){
read_u(u0, K1, K2, initfile);
return 0;
}

5
src/init.h
View File

@@ -1,10 +1,15 @@
#ifndef INIT_H
#define INIT_H
#include <stdio.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);
// Initialize from file
int init_file (_Complex double* u0, int K1, int K2, FILE* initfile);
#endif

112
src/io.c Normal file
View File

@@ -0,0 +1,112 @@
#include <stdbool.h>
#include <stdlib.h>
#include "io.h"
#include "navier-stokes.h"
// write final entry to file
int write_u(_Complex double* u, int K1, int K2, FILE* file){
int kx,ky;
// do nothing if there is no file
if(file==NULL){
return 0;
}
for(kx=-K1;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){
fprintf(file,"%d %d % .15e % .15e\n",kx,ky,__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)],__imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]);
}
}
return 0;
}
// read u from file
int read_u(_Complex double* u, int K1, int K2, FILE* file){
int kx,ky;
double r,i;
char* line;
unsigned int len=256;
unsigned int pos=0;
char* line_realloc;
char c;
int ret;
unsigned int counter=0;
bool commented=false;
// error if there is no file (this should not happen)
if (file==NULL){
fprintf(stderr,"error reading u from file (this is a bug!)\n");
return -1;
}
// allocate line buffer
line=calloc(sizeof(char), len);
while(1){
c=fgetc(file);
// end of file
if (feof(file)){
break;
}
// newline: read line and reset buffer
if(c=='\n' || c=='\r'){
// increment line counter
counter++;
// read entry
// ignore empty lines
if(pos>0){
ret=sscanf(line, "%d %d %le %le", &kx, &ky, &r, &i);
// errors
if(ret!=4){
fprintf(stderr, "warning: line %d does not match the input format: '%s'\n", counter, line);
}
else{
if(kx>K1 || kx<-K1 || ky>K2 || ky<-K2){
fprintf(stderr, "warning: reading line %d: kx or ky out of bounds: %d,%d\n", counter, kx, ky);
}
else{
// set u
u[klookup(kx, ky, 2*K1+1, 2*K2+1)]=r+i*I;
}
}
}
// reset buffer
pos=0;
line[pos]='\0';
commented=false;
}
// comment: stop reading
else if (c=='#'){
commented=true;
}
// add to buffer (unless we are in a comment)
else if (!(commented)){
// check that there is room in buffer
if(pos==len){
// too short: reallocate
line_realloc=calloc(sizeof(char), 2*len);
for(pos=0;pos<len;pos++){
line_realloc[pos]=line[pos];
}
free(line);
line=line_realloc;
len=2*len;
}
// add c to line
line[pos]=c;
pos++;
line[pos]='\0';
}
}
free(line);
return 0;
}

12
src/io.h Normal file
View File

@@ -0,0 +1,12 @@
#ifndef IO_H
#define IO_H
#include <stdio.h>
// write u to file
int write_u(_Complex double* u, int K1, int K2, FILE* file);
// read u from file
int read_u(_Complex double* u, int K1, int K2, FILE* file);
#endif

83
src/main.c
View File

@@ -6,6 +6,7 @@
#include <string.h>
#include <stdlib.h>
#include <stdbool.h>
#include <errno.h>
#include "navier-stokes.h"
#include "driving.h"
#include "init.h"
@@ -28,17 +29,17 @@ typedef struct nstrophy_parameters {
// usage message
int print_usage();
// print parameters
int print_params(nstrophy_parameters parameters, unsigned int driving, unsigned int init, FILE* file);
int print_params(nstrophy_parameters parameters, unsigned int driving, unsigned int init, char* initfile_str, FILE* file);
// read command line arguments
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_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* init, unsigned int* nthreads, char** savefile_str, char** initfile_str);
int read_params(char* param_str, nstrophy_parameters* parameters);
int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2);
// set driving force
_Complex double* set_driving(unsigned int driving, nstrophy_parameters parameters);
// set initial condition
_Complex double* set_init(unsigned int init, nstrophy_parameters parameters);
_Complex double* set_init(unsigned int init, nstrophy_parameters parameters, FILE* initfile);
#define COMMAND_UK 1
#define COMMAND_EEA 2
@@ -49,6 +50,7 @@ _Complex double* set_init(unsigned int init, nstrophy_parameters parameters);
#define INIT_RANDOM 1
#define INIT_GAUSSIAN 2
#define INIT_FILE 3
int main (
@@ -62,13 +64,17 @@ int main (
unsigned int nthreads=1;
_Complex double* u0;
_Complex double *g;
char* savefile_str=NULL;
char* initfile_str=NULL;
FILE* savefile=NULL;
FILE* initfile=NULL;
command=0;
driving=0;
init=0;
// read command line arguments
ret=read_args(argc, argv, &param_str, &driving, &init, &command, &nthreads);
ret=read_args(argc, argv, &param_str, &driving, &init, &command, &nthreads, &savefile_str, &initfile_str);
if(ret<0){
return(-1);
}
@@ -79,24 +85,47 @@ int main (
return(-1);
}
// open initfile
if(initfile_str!=NULL){
initfile=fopen(initfile_str,"r");
if(initfile==NULL){
fprintf(stderr,"Error opening file '%s' for reading: %s\n", initfile_str, strerror(errno));
return(-1);
}
}
// set driving force
g=set_driving(driving, parameters);
// set initial condition
u0=set_init(init, parameters);
u0=set_init(init, parameters, initfile);
// close initfile (do this early, so that it is possible to use the same file for init and save)
if (initfile!=NULL){
fclose(initfile);
}
// open savefile (do this after closing init file)
if(savefile_str!=NULL){
savefile=fopen(savefile_str,"w");
if(savefile==NULL){
fprintf(stderr,"Error opening file '%s' for writing: %s\n", savefile_str, strerror(errno));
return(-1);
}
}
// print parameters
print_params(parameters, driving, init, stderr);
print_params(parameters, driving, init, stdout);
print_params(parameters, driving, init, initfile_str, stderr);
print_params(parameters, driving, init, initfile_str, stdout);
// run command
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, nthreads);
uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads, savefile);
}
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, nthreads);
eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads, savefile);
}
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);
quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, nthreads, savefile);
}
else if(command==0){
fprintf(stderr, "error: no command specified\n");
@@ -106,6 +135,11 @@ int main (
free(g);
free(u0);
// close savefile
if (savefile!=NULL){
fclose(savefile);
}
return(0);
}
@@ -120,6 +154,7 @@ int print_params(
nstrophy_parameters parameters,
unsigned int driving,
unsigned int init,
char* initfile_str,
FILE* file
){
fprintf(file,"# K1=%d, K2=%d, N1=%d, N2=%d, nu=%.15e, delta=%.15e, L=%.15e", parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nu, parameters.delta, parameters.L);
@@ -143,6 +178,9 @@ int print_params(
case INIT_GAUSSIAN:
fprintf(file,", init=gaussian");
break;
case INIT_FILE:
fprintf(file,", init=file:%s", initfile_str);
break;
default:
fprintf(file,", init=gaussian");
break;
@@ -158,6 +196,7 @@ int print_params(
#define CP_FLAG_DRIVING 2
#define CP_FLAG_NTHREADS 3
#define CP_FLAG_INIT 4
#define CP_FLAG_SAVEFILE 5
int read_args(
int argc,
const char* argv[],
@@ -165,7 +204,9 @@ int read_args(
unsigned int* driving_force,
unsigned int* init,
unsigned int* command,
unsigned int* nthreads
unsigned int* nthreads,
char** savefile_str,
char** initfile_str
){
int i;
int ret;
@@ -192,6 +233,9 @@ int read_args(
case 'i':
flag=CP_FLAG_INIT;
break;
case 's':
flag=CP_FLAG_SAVEFILE;
break;
default:
fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
print_usage();
@@ -236,12 +280,22 @@ int read_args(
else if (strcmp(argv[i],"gaussian")==0){
*init=INIT_GAUSSIAN;
}
// matches any argument that starts with 'file:'
else if (strncmp(argv[i],"file:",5)==0){
*init=INIT_FILE;
*initfile_str=(char*)argv[i]+5;
}
else{
fprintf(stderr, "error: unrecognized initial condition '%s'\n",argv[i]);
return(-1);
}
flag=0;
}
// savefile
else if(flag==CP_FLAG_SAVEFILE){
*savefile_str=(char*)argv[i];
flag=0;
}
// computation to run
else{
if(strcmp(argv[i], "uk")==0){
@@ -495,7 +549,8 @@ _Complex double* set_driving(
// set initial condition
_Complex double* set_init(
unsigned int init,
nstrophy_parameters parameters
nstrophy_parameters parameters,
FILE* initfile
){
_Complex double* u0=calloc(sizeof(_Complex double),(2*parameters.K1+1)*(2*parameters.K2+1));
@@ -508,6 +563,10 @@ _Complex double* set_init(
init_gaussian(u0, parameters.K1, parameters.K2);
break;
case INIT_FILE:
init_file(u0, parameters.K1, parameters.K2, initfile);
break;
default:
init_gaussian(u0, parameters.K1, parameters.K2);
break;

19
src/navier-stokes.c
View File

@@ -1,4 +1,5 @@
#include "navier-stokes.h"
#include "io.h"
#include <math.h>
#include <stdlib.h>
@@ -15,7 +16,8 @@ int uk(
_Complex double* u0,
_Complex double* g,
unsigned int print_freq,
unsigned int nthreads
unsigned int nthreads,
FILE* savefile
){
_Complex double* u;
_Complex double* tmp1;
@@ -65,6 +67,9 @@ int uk(
}
}
// save final entry to savefile
write_u(u, K1, K2, savefile);
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
return(0);
}
@@ -82,7 +87,8 @@ int eea(
_Complex double* u0,
_Complex double* g,
unsigned int print_freq,
unsigned int nthreads
unsigned int nthreads,
FILE* savefile
){
_Complex double* u;
_Complex double* tmp1;
@@ -126,6 +132,9 @@ int eea(
}
}
// save final entry to savefile
write_u(u, K1, K2, savefile);
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
return(0);
}
@@ -142,7 +151,8 @@ int quiet(
double L,
_Complex double* u0,
_Complex double* g,
unsigned int nthreads
unsigned int nthreads,
FILE* savefile
){
_Complex double* u;
_Complex double* tmp1;
@@ -162,6 +172,9 @@ int quiet(
ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
}
// save final entry to savefile
write_u(u, K1, K2, savefile);
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
return(0);
}

9
src/navier-stokes.h
View File

@@ -13,16 +13,13 @@ typedef struct fft_vects {
} fft_vect;
// 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 nthreads);
// 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* u0, _Complex double* g, 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, unsigned int print_freq, unsigned int nthreadsl, FILE* savefile);
// 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 nthreads);
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 nthreads, FILE* savefile);
// 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);
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);
// initialize vectors for computation