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

write ns_step routine that handles all the algorithms

Browse Source
This commit is contained in:
Ian Jauslin
2024-02-19 14:42:34 -05:00
parent 9fe9e3d96f
commit 9ecf4413a5
2 changed files with 63 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files

108
src/navier-stokes.c
View File

@@ -18,7 +18,6 @@ limitations under the License.
#include "io.h"
#include "navier-stokes.h"
#include "complex_tools.h"
#include <cblas.h>
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
@@ -90,21 +89,8 @@ int uk(
// iterate
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RKF45) {
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
} else if (algorithm==ALGORITHM_RKDP54) {
// only compute k1 on the first step
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
} else if (algorithm==ALGORITHM_RKBS32) {
// only compute k1 on the first step
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
// step
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
time+=step;
step=next_step;
@@ -206,21 +192,7 @@ int enstrophy(
// iterate
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RKF45) {
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
} else if (algorithm==ALGORITHM_RKDP54) {
// only compute k1 on the first step
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
} else if (algorithm==ALGORITHM_RKBS32) {
// only compute k1 on the first step
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
time+=step;
@@ -387,21 +359,7 @@ int quiet(
// iterate
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RKF45) {
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
} else if (algorithm==ALGORITHM_RKDP54) {
// only compute k1 on the first step
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
} else if (algorithm==ALGORITHM_RKBS32) {
// only compute k1 on the first step
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
time+=step;
step=next_step;
@@ -554,6 +512,60 @@ int copy_u(
}
// next time step
int ns_step(
unsigned int algorithm,
_Complex double* u,
int K1,
int K2,
int N1,
int N2,
double nu,
double* delta,
double* next_delta,
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
double L,
_Complex double* g,
double time,
double starting_time,
fft_vect fft1,
fft_vect fft2,
fft_vect ifft,
// the pointers tmp1 and tmp2 will be exchanged at the end of the routine for first-same-as-last RK algorithms
_Complex double** tmp1,
_Complex double** tmp2,
_Complex double* tmp3,
_Complex double* tmp4,
_Complex double* tmp5,
_Complex double* tmp6,
_Complex double* tmp7,
bool irreversible,
bool keep_en_cst,
double target_en
){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, *delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, irreversible, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, *delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, tmp3, irreversible, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RKF45) {
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
} else if (algorithm==ALGORITHM_RKDP54) {
// only compute k1 on the first step
// first-same-as-last with 2-nd argument
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
} else if (algorithm==ALGORITHM_RKBS32) {
// only compute k1 on the first step
// first-same-as-last with 4-th argument
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp3, tmp4, tmp2, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers.edu\n",algorithm);
}
return (0);
}
// RK 4 algorithm
int ns_step_rk4(
_Complex double* u,
@@ -1343,7 +1355,7 @@ int ns_T_nofft(
return 0;
}
/*
// Jacobian of rhs
int ns_D_rhs(
_Complex double* out,
@@ -1365,7 +1377,6 @@ int ns_D_rhs(
alpha=0;
}
#define MATSIZE (K1*(2*(K2+1)+K2))
for(i=0; i<MATSIZE*MATSIZE; i++){
out[i]=0;
}
@@ -1414,6 +1425,7 @@ _Complex double ns_d_T(
return (ky*lx-kx*ly)*(qx*qx+qy*qy-lx*lx-ly*ly)/sqrt(lx*lx+ly*ly)/sqrt(qx*qx+qy*qy)*u[klookup_sym(-qx,-qy,K2)];
}
}
*/
// compute alpha