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

Driving force as vector instead of function

Browse Source
This commit is contained in:
Ian Jauslin 2022-05-26 15:16:44 -04:00
parent d4254c6b8e
commit aa66aadb74
5 changed files with 76 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/driving.c
View File

@ -1,28 +1,40 @@
#include "driving.h"
#include "navier-stokes.h"
#include <math.h>
_Complex double g_test(
int kx,
int ky
// test driving function
int g_test(
_Complex double* g,
int K1,
int K2
){
//return sqrt(kx*kx*ky*ky)*exp(-(kx*kx+ky*ky));
int kx,ky;
for(kx=-K1;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){
if(kx==2 && ky==-1){
return 0.5+sqrt(3)/2*I;
g[klookup(kx,ky,K1,K2)]=0.5+sqrt(3)/2*I;
}
else if(kx==-2 && ky==1){
return 0.5-sqrt(3)/2*I;
g[klookup(kx,ky,K1,K2)]=0.5-sqrt(3)/2*I;
}
else{
g[klookup(kx,ky,K1,K2)]=0.;
}
}
return 0.;
}
#define UNUSED(x) (void)(x)
_Complex double g_zero(
int kx,
int ky
return 0;
}
int g_zero(
_Complex double* g,
int K1,
int K2
){
// avoid unused variable warnings
UNUSED(kx);
UNUSED(ky);
return 0.;
int i;
for(i=0;i<(2*K1+1)*(2*K2+1);i++){
g[i]=0.;
}
return 0;
}

6
src/driving.h
View File

@ -1,9 +1,7 @@
#ifndef DRIVING_H
#define DRIVING_H
#include <complex.h>
_Complex double g_zero( int kx, int ky);
_Complex double g_test( int kx, int ky);
int g_zero(_Complex double* g, int K1, int K2);
int g_test(_Complex double* g, int K1, int K2);
#endif

41
src/main.c
View File

@ -32,6 +32,8 @@ int read_args(int argc, const char* argv[], char** params, unsigned int* driving
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);
@ -53,11 +55,11 @@ int main (
){
char* param_str=NULL;
nstrophy_parameters parameters;
_Complex double (*g)(int,int);
int ret;
unsigned int driving,command,init;
unsigned int nthreads=1;
_Complex double* u0;
_Complex double *g;
command=0;
driving=0;
@ -76,19 +78,7 @@ int main (
}
// set driving force
switch(driving){
case DRIVING_ZERO:
g=g_zero;
break;
case DRIVING_TEST:
g=g_test;
break;
default:
g=g_zero;
break;
}
g=set_driving(driving, parameters);
// set initial condition
u0=set_init(init, parameters);
@ -110,6 +100,7 @@ int main (
print_usage();
}
free(g);
free(u0);
return(0);
@ -441,6 +432,28 @@ int set_parameter(
return(0);
}
// set driving force
_Complex double* set_driving(
unsigned int driving,
nstrophy_parameters parameters
){
_Complex double* g=calloc(sizeof(_Complex double),(2*parameters.K1+1)*(2*parameters.K2+1));
switch(driving){
case DRIVING_ZERO:
g_zero(g, parameters.K1, parameters.K2);
break;
case DRIVING_TEST:
g_test(g, parameters.K1, parameters.K2);
break;
default:
g_test(g, parameters.K1, parameters.K2);
break;
}
return g;
}
// set initial condition
_Complex double* set_init(
unsigned int init,

18
src/navier-stokes.c
View File

@ -13,7 +13,7 @@ int uk(
double delta,
double L,
_Complex double* u0,
_Complex double (*g)(int,int),
_Complex double* g,
unsigned int print_freq,
unsigned int nthreads
){
@ -80,7 +80,7 @@ int energy(
double delta,
double L,
_Complex double* u0,
_Complex double (*g)(int,int),
_Complex double* g,
unsigned int print_freq,
unsigned int nthreads
){
@ -131,7 +131,7 @@ int enstrophy(
double delta,
double L,
_Complex double* u0,
_Complex double (*g)(int,int),
_Complex double* g,
unsigned int print_freq,
unsigned int nthreads
){
@ -185,7 +185,7 @@ int quiet(
double delta,
double L,
_Complex double* u0,
_Complex double (*g)(int,int),
_Complex double* g,
unsigned int nthreads
){
_Complex double* u;
@ -306,7 +306,7 @@ int ins_step(
double nu,
double delta,
double L,
_Complex double (*g)(int,int),
_Complex double* g,
fft_vect fft1,
fft_vect fft2,
fft_vect ifft,
@ -383,7 +383,7 @@ int ins_rhs(
int N2,
double nu,
double L,
_Complex double (*g)(int,int),
_Complex double* g,
fft_vect fft1,
fft_vect fft2,
fft_vect ifft
@ -413,7 +413,7 @@ int ins_rhs(
for(ky=-K2;ky<=K2;ky++){
if(kx!=0 || ky!=0){
// enforce the reality of u by adding ifft.fft(k) and the conjugate of ifft.fft(-k)
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)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
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[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
}
}
}
@ -535,7 +535,7 @@ _Complex double compute_alpha(
_Complex double* u,
int K1,
int K2,
_Complex double (*g)(int,int)
_Complex double* g
){
_Complex double num=0;
_Complex double denom=0;
@ -544,7 +544,7 @@ _Complex double compute_alpha(
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)])*(1+(ky!=0?kx*kx/ky/ky:0));
num+=(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj((*g)(kx,ky))*(1+(ky!=0?kx*kx/ky/ky:0));
num+=(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(g[klookup(kx,ky,2*K1+1,2*K2+1)])*(1+(ky!=0?kx*kx/ky/ky:0));
}
}

14
src/navier-stokes.h
View File

@ -13,16 +13,16 @@ 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)(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, 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)(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, unsigned int print_freq, unsigned int nthreads);
// compute enstrophy
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);
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, unsigned int print_freq, unsigned int nthreads);
// 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)(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, unsigned int nthreads);
// initialize vectors for computation
@ -34,10 +34,10 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);
// 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, 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
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);
int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft);
// 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);
@ -46,7 +46,7 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
// compute alpha
_Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double (*g)(int,int));
_Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double* g);
// get index for kx,ky in array of size S
int klookup( int kx, int ky, int S1, int S2);