diff options
| author | Ian Jauslin <ian.jauslin@rutgers.edu> | 2023-04-26 11:13:50 -0400 |
|---|---|---|
| committer | Ian Jauslin <ian.jauslin@rutgers.edu> | 2023-04-26 11:13:50 -0400 |
| commit | 527365d62e8c770142c3c1065b5f973433bd60b2 (patch) | |
| tree | 5adabaa1f96c20be42dab7b0e9602c7e0177945c /src/navier-stokes.c | |
| parent | f7a7a5866ce3c43f497e13734d968b9292254014 (diff) | |
Implement RK2
Diffstat (limited to 'src/navier-stokes.c')
| -rw-r--r-- | src/navier-stokes.c | 50 |
1 files changed, 46 insertions, 4 deletions
diff --git a/src/navier-stokes.c b/src/navier-stokes.c index ce07077..755e829 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -51,7 +51,7 @@ int uk( // iterate 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){ fprintf(stderr,"%lu % .8e ",t,t*delta); @@ -127,7 +127,7 @@ int eea( // iterate 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); alpha=compute_alpha(u, K1, K2, g, L); @@ -214,7 +214,7 @@ int quiet( // iterate 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 @@ -311,7 +311,8 @@ int copy_u( } // next time step -int ns_step( +// RK 4 algorithm +int ns_step_rk4( _Complex double* u, int K1, int K2, @@ -388,6 +389,47 @@ int ns_step( 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 int ns_rhs( _Complex double* out, |