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

Fix initial enstrophy

Browse Source
This commit is contained in:
Ian Jauslin
2023-04-05 22:41:19 -04:00
parent 0bf223bcb9
commit 1616b6bbae
3 changed files with 53 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
src/init.c
View File

@@ -7,9 +7,12 @@
// random initial condition
int init_random (
_Complex double* u0,
double init_en,
int K1,
int K2,
int seed
double L,
int seed,
bool irreversible
){
int kx,ky;
double rescale;
@@ -29,16 +32,16 @@ int init_random (
}
}
// 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);
}
if (irreversible) {
// fix the energy
rescale=compute_energy(u0, K1, K2);
} else {
// fix the enstrophy
rescale=compute_enstrophy(u0, K1, K2, L);
}
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);
u0[klookup(kx,ky,2*K1+1,2*K2+1)]=u0[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(init_en/rescale);
}
}
@@ -48,10 +51,14 @@ int init_random (
// Gaussian initial condition
int init_gaussian (
_Complex double* u0,
double init_en,
int K1,
int K2
int K2,
double L,
bool irreversible
){
int kx,ky;
double rescale;
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
@@ -59,6 +66,19 @@ int init_gaussian (
}
}
if (irreversible) {
// fix the energy
rescale=compute_energy(u0, K1, K2);
} else {
// fix the enstrophy
rescale=compute_enstrophy(u0, K1, K2, L);
}
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(init_en/rescale);
}
}
return 0;
}

5
src/init.h
View File

@@ -2,12 +2,13 @@
#define INIT_H
#include <stdio.h>
#include <stdbool.h>
// random initial condition
int init_random(_Complex double* u0, int K1, int K2, int seed);
int init_random(_Complex double* u0, double init_en, int K1, int K2, double L, int seed, bool irreversible);
// Gaussian initial condition
int init_gaussian(_Complex double* u0, int K1, int K2);
int init_gaussian(_Complex double* u0, double init_en, int K1, int K2, double L, bool irreversible);
// Initialize from file
int init_file (_Complex double* u0, int K1, int K2, FILE* initfile);

25
src/main.c
View File

@@ -14,6 +14,7 @@
// structure to store parameters, to make it easier to pass parameters to CLI functions
typedef struct nstrophy_parameters {
double init_en; // initial value for the energy for ins and enstrophy for rns
bool irreversible;
int K1;
int K2;
@@ -159,7 +160,15 @@ int print_params(
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);
fprintf(file,"# ");
if (parameters.irreversible){
fprintf(file,"equation=irreversible");
} else {
fprintf(file,"equation=reversible");
}
fprintf(file,", K1=%d, K2=%d, N1=%d, N2=%d, nu=%.15e, delta=%.15e, L=%.15e, init_en=%.15e", parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nu, parameters.delta, parameters.L, parameters.init_en);
switch(driving){
case DRIVING_TEST:
@@ -338,6 +347,7 @@ int read_params(
bool lhs=true;
// defaults
parameters->init_en=1.54511597324389e+02;
parameters->irreversible=true;
parameters->K1=16;
parameters->K2=parameters->K1;
@@ -442,6 +452,13 @@ int set_parameter(
return(-1);
}
}
else if (strcmp(lhs,"init_en")==0){
ret=sscanf(rhs,"%lf",&(parameters->init_en));
if(ret!=1){
fprintf(stderr, "error: parameter 'init_en' should be a double\n got '%s'\n",rhs);
return(-1);
}
}
else if (strcmp(lhs,"K1")==0){
ret=sscanf(rhs,"%d",&(parameters->K1));
if(ret!=1){
@@ -579,11 +596,11 @@ _Complex double* set_init(
switch(init){
case INIT_RANDOM:
init_random(u0, parameters.K1, parameters.K2, parameters.seed);
init_random(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.seed, parameters.irreversible);
break;
case INIT_GAUSSIAN:
init_gaussian(u0, parameters.K1, parameters.K2);
init_gaussian(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.irreversible);
break;
case INIT_FILE:
@@ -591,7 +608,7 @@ _Complex double* set_init(
break;
default:
init_gaussian(u0, parameters.K1, parameters.K2);
init_gaussian(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.irreversible);
break;
}