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

Fix running averages for adaptive step, and base print_freq and final_time on times

Browse Source
This commit is contained in:
Ian Jauslin
2023-05-16 00:00:18 -04:00
parent 5db3680b31
commit 5ede4f1084
6 changed files with 90 additions and 157 deletions
Download Patch File Download Diff File Expand all files Collapse all files

109
src/navier-stokes.c
View File

@@ -17,8 +17,8 @@ limitations under the License.
#include "constants.cpp"
#include "io.h"
#include "navier-stokes.h"
#include "statistics.h"
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
@@ -28,7 +28,7 @@ int uk(
int K2,
int N1,
int N2,
uint64_t nsteps,
double final_time,
double nu,
double delta,
double L,
@@ -38,13 +38,11 @@ int uk(
_Complex double* g,
bool irreversible,
unsigned int algorithm,
uint64_t print_freq,
uint64_t starting_step,
double print_freq,
double starting_time,
unsigned int nthreads,
FILE* savefile
){
double time=starting_time;
_Complex double* u;
_Complex double* tmp1;
_Complex double* tmp2;
@@ -53,7 +51,7 @@ int uk(
_Complex double* tmp5;
_Complex double* tmp6;
_Complex double* tmp7;
uint64_t t;
double time;
fft_vect fft1;
fft_vect fft2;
fft_vect ifft;
@@ -65,18 +63,23 @@ int uk(
// print column headers
printf("# 1:i 2:t ");
t=3;
unsigned int i=3;
for(kx=-K1;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){
printf(" %6lu:(%4d,%4d)r ",t,kx,ky);
t++;
printf(" %6lu:(%4d,%4d)i ",t,kx,ky);
t++;
printf(" %6u:(%4d,%4d)r ",i,kx,ky);
i++;
printf(" %6u:(%4d,%4d)i ",i,kx,ky);
i++;
}
}
// period
// add 0.1 to ensure proper rounding
uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1);
// iterate
for(t=starting_step;nsteps==0 || t<starting_step+nsteps;t++){
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else if (algorithm==ALGORITHM_RK4) {
@@ -89,9 +92,11 @@ int uk(
time+=delta;
if(t%print_freq==0){
fprintf(stderr,"%lu % .8e ",t,time);
printf("%8lu % .15e ",t,time);
if(time>(n+1)*print_freq){
n++;
fprintf(stderr,"% .8e ",time);
printf("% .15e ",time);
for(kx=-K1;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){
@@ -120,7 +125,7 @@ int enstrophy(
int K2,
int N1,
int N2,
uint64_t nsteps,
double final_time,
double nu,
double delta,
double L,
@@ -130,8 +135,7 @@ int enstrophy(
_Complex double* g,
bool irreversible,
unsigned int algorithm,
uint64_t print_freq,
uint64_t starting_step,
double print_freq,
double starting_time,
unsigned int nthreads,
FILE* savefile,
@@ -148,11 +152,11 @@ int enstrophy(
_Complex double* tmp5;
_Complex double* tmp6;
_Complex double* tmp7;
double time=starting_time;
double time;
double alpha, enstrophy;
double prevtime;
double avg_a,avg_en,avg_en_x_a;
// index
uint64_t t;
fft_vect fft1;
fft_vect fft2;
fft_vect ifft;
@@ -166,12 +170,15 @@ int enstrophy(
avg_a=0;
avg_en=0;
avg_en_x_a=0;
prevtime=starting_time;
// special first case when starting_time is not a multiple of print_freq
uint64_t first_box = print_freq - (starting_step % print_freq);
// period
// add 0.1 to ensure proper rounding
uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1);
// iterate
for(t=starting_step;nsteps==0 || t<starting_step+nsteps;t++){
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else if (algorithm==ALGORITHM_RK4) {
@@ -187,21 +194,35 @@ int enstrophy(
alpha=compute_alpha(u, K1, K2, g, L);
enstrophy=compute_enstrophy(u, K1, K2, L);
avg_a=average_step(alpha, avg_a, t, starting_step, print_freq, first_box);
avg_en=average_step(enstrophy, avg_en, t, starting_step, print_freq, first_box);
avg_en_x_a=average_step(enstrophy*alpha, avg_en_x_a, t, starting_step, print_freq, first_box);
// add to running averages (estimate the total duration of interval as print_freq, will be adjusted later)
avg_a+=alpha*(delta/print_freq);
avg_en+=enstrophy*(delta/print_freq);
avg_en_x_a+=enstrophy*alpha*(delta/print_freq);
if(time>(n+1)*print_freq){
n++;
// adjust duration of interval to its actual value
avg_a*=print_freq/(time-prevtime);
avg_en*=print_freq/(time-prevtime);
avg_en_x_a*=print_freq/(time-prevtime);
if(t>starting_step && t%print_freq==0){
// print to stderr so user can follow along
if(algorithm==ALGORITHM_RKF45){
fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, delta);
fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, delta);
} else {
fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
}
// print to stdout
printf("%8lu % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
}
// reset averages
avg_a=0;
avg_en=0;
avg_en_x_a=0;
prevtime=time;
// catch abort signal
if (g_abort){
// print u to stderr if no savefile
@@ -223,22 +244,20 @@ int enstrophy(
if(params_string!=NULL) {
char* params=calloc(sizeof(char), strlen(params_string)+1);
strcpy(params, params_string);
remove_entry(params, "starting_step");
remove_entry(params, "starting_time");
remove_entry(params, "init");
remove_entry(params, "nsteps");
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
remove_entry(params, "delta");
}
fprintf(savefile," -p \"%s;starting_step=%lu;nsteps=%lu;init=file:%s", params, t+1, (nsteps+starting_step < t+1 ? 0 : nsteps+starting_step-t-1), savefile_string);
fprintf(savefile," -p \"%s;init=file:%s", params, savefile_string);
free(params);
// write delta if adaptive, and not writing binary
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){
fprintf(savefile,";delta=%.15e;starting_time=%.15e", delta, starting_time);
fprintf(savefile,";delta=%.15e", delta);
}
// write starting_time if not adaptive
if(algorithm<ALGORITHM_ADAPTIVE_THRESHOLD){
fprintf(savefile,";starting_time=%.15e", starting_time);
// write starting_time if not writing binary
if(savefile==stderr || savefile==stdout){
fprintf(savefile,";starting_time=%.15e", time);
}
fprintf(savefile,"\"");
}
@@ -249,12 +268,11 @@ int enstrophy(
write_vec(u, K1, K2, savefile);
} else {
write_vec_bin(u, K1, K2, savefile);
// last binary entry: starting time
fwrite(&time, sizeof(double), 1, savefile);
// extra binary data for adaptive algorithm
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
// first binary entry: delta
fwrite(&delta, sizeof(double), 1, savefile);
// ssecond binary entry: starting time
fwrite(&starting_time, sizeof(double), 1, savefile);
}
}
}
@@ -269,13 +287,13 @@ int quiet(
int K2,
int N1,
int N2,
uint64_t nsteps,
double final_time,
double nu,
double delta,
double L,
double adaptive_tolerance,
double adaptive_factor,
uint64_t starting_step,
double starting_time,
_Complex double* u0,
_Complex double* g,
bool irreversible,
@@ -291,7 +309,7 @@ int quiet(
_Complex double* tmp5;
_Complex double* tmp6;
_Complex double* tmp7;
uint64_t t;
double time;
fft_vect fft1;
fft_vect fft2;
fft_vect ifft;
@@ -301,7 +319,8 @@ int quiet(
copy_u(u, u0, K1, K2);
// iterate
for(t=starting_step;nsteps==0 || t<starting_step+nsteps;t++){
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else if (algorithm==ALGORITHM_RK4) {
@@ -311,6 +330,8 @@ int quiet(
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
time+=delta;
}
// save final entry to savefile