Implement RK2

This commit is contained in:
Ian Jauslin 2023-04-26 11:13:50 -04:00
parent f7a7a5866c
commit 527365d62e
2 changed files with 50 additions and 5 deletions

View File

@ -51,7 +51,7 @@ int uk(
// iterate // iterate
for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){ for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
if(t%print_freq==0){ if(t%print_freq==0){
fprintf(stderr,"%lu % .8e ",t,t*delta); fprintf(stderr,"%lu % .8e ",t,t*delta);
@ -127,7 +127,7 @@ int eea(
// iterate // iterate
for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){ for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
energy=compute_energy(u, K1, K2); energy=compute_energy(u, K1, K2);
alpha=compute_alpha(u, K1, K2, g, L); alpha=compute_alpha(u, K1, K2, g, L);
@ -214,7 +214,7 @@ int quiet(
// iterate // iterate
for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){ for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
} }
// save final entry to savefile // save final entry to savefile
@ -311,7 +311,8 @@ int copy_u(
} }
// next time step // next time step
int ns_step( // RK 4 algorithm
int ns_step_rk4(
_Complex double* u, _Complex double* u,
int K1, int K1,
int K2, int K2,
@ -388,6 +389,47 @@ int ns_step(
return(0); return(0);
} }
// RK 2 algorithm
int ns_step_rk2(
_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,
bool irreversible
){
int kx,ky;
// k1
ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// u+h*k1/2
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
}
}
// k2
ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
u[klookup_sym(kx,ky,K2)]+=delta*tmp1[klookup_sym(kx,ky,K2)];
}
}
return(0);
}
// right side of Irreversible/Reversible Navier-Stokes equation // right side of Irreversible/Reversible Navier-Stokes equation
int ns_rhs( int ns_rhs(
_Complex double* out, _Complex double* out,

View File

@ -36,7 +36,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); int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);
// next time step for Irreversible/reversible Navier-Stokes equation // next time step for Irreversible/reversible Navier-Stokes equation
int ns_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, bool irreversible); // RK4
int ns_step_rk4( _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, bool irreversible);
// RK2
int ns_step_rk2( _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, bool irreversible);
// right side of Irreversible/reversible Navier-Stokes equation // right side of Irreversible/reversible Navier-Stokes equation
int ns_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, bool irreversible); int ns_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, bool irreversible);