Jacobian of un->un+1

This commit is contained in:
Ian Jauslin 2024-02-19 15:20:19 -05:00
parent 9ecf4413a5
commit 89601d19d1
5 changed files with 134 additions and 2 deletions

View File

@ -14,6 +14,7 @@ See the License for the specific language governing permissions and
limitations under the License. limitations under the License.
*/ */
#define M_PI 3.14159265358979323846
#define COMMAND_UK 1 #define COMMAND_UK 1
#define COMMAND_ENSTROPHY 2 #define COMMAND_ENSTROPHY 2

107
src/lyapunov.c Normal file
View File

@ -0,0 +1,107 @@
/*
Copyright 2017-2023 Ian Jauslin
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
#include "constants.cpp"
#include "lyapunov.h"
#include <cblas.h>
#define MATSIZE (K1*(2*(K2+1)+K2))
// Jacobian of u_n -> u_{n+1}
int ns_D_step(
_Complex double* out,
_Complex double* un,
_Complex double* un_next,
int K1,
int K2,
int N1,
int N2,
double nu,
double epsilon,
double delta,
unsigned int algorithm,
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
double L,
_Complex double* g,
double time,
fft_vect fft1,
fft_vect fft2,
fft_vect ifft,
_Complex double* tmp1,
_Complex double* tmp2,
_Complex double* tmp3,
_Complex double* tmp4,
_Complex double* tmp5,
_Complex double* tmp6,
_Complex double* tmp7,
_Complex double* tmp8,
bool irreversible,
bool keep_en_cst,
double target_en
){
int lx,ly;
int i;
double step, next_step;
for(lx=0;lx<=K1;lx++){
for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
// derivative in the real direction
// finite difference vector
for (i=0;i<MATSIZE;i++){
if(i==klookup_sym(lx,ly,K2)){
tmp1[i]=un[i]+epsilon;
}else{
tmp1[i]=un[i];
}
}
// compute step
// reinitialize step
step=delta;
next_step=delta;
ns_step(algorithm, tmp1, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, time, fft1, fft2, ifft, &tmp2, &tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, irreversible, keep_en_cst, target_en);
// compute derivative
for (i=0;i<MATSIZE;i++){
out[i]=(tmp1[i]-un_next[i])/epsilon;
}
// derivative in the imaginary direction
// finite difference vector
for (i=0;i<MATSIZE;i++){
if(i==klookup_sym(lx,ly,K2)){
tmp1[i]=un[i]+epsilon*I;
}else{
tmp1[i]=un[i];
}
}
// compute step
// reinitialize step
step=delta;
next_step=delta;
ns_step(algorithm, tmp1, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, time, fft1, fft2, ifft, &tmp2, &tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, irreversible, keep_en_cst, target_en);
// compute derivative
for (i=0;i<MATSIZE;i++){
out[i]+=-(tmp1[i]-un_next[i])/epsilon*I;
}
}
}
return(0);
}

25
src/lyapunov.h Normal file
View File

@ -0,0 +1,25 @@
/*
Copyright 2017-2023 Ian Jauslin
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
#ifndef LYAPUNOV_H
#define LYAPUNOV_H
#include "navier-stokes.h"
int ns_D_step( _Complex double* out, _Complex double* un, _Complex double* un_next, int K1, int K2, int N1, int N2, double nu, double epsilon, double delta, unsigned int algorithm, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double L, _Complex double* g, double time, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, _Complex double* tmp8, bool irreversible, bool keep_en_cst, double target_en);
#endif

View File

@ -1355,6 +1355,7 @@ int ns_T_nofft(
return 0; return 0;
} }
/* /*
// Jacobian of rhs // Jacobian of rhs
int ns_D_rhs( int ns_D_rhs(

View File

@ -21,8 +21,6 @@ limitations under the License.
#include <stdbool.h> #include <stdbool.h>
#include <fftw3.h> #include <fftw3.h>
#define M_PI 3.14159265358979323846
// abort signal // abort signal
extern volatile bool g_abort; extern volatile bool g_abort;