我试图将一个函数
slice
的输出切片传递给另一个函数,但是我在 append_2d_slice
切片函数的第 93 行收到分段错误错误。我评论了导致问题的行。我怀疑我以错误的方式将 slice
变量传递给函数,但我不确定。
这是代码:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//Works on Collab
void awe_2d_heterogeneous_8th_order_data_time_reverse(
double** UUo, double** UUm, int nx, int ny, double dx, double dy,
double dt, double** v, double** slice, int it, double rx
) {
// Get dimensions of the wavefield
int irx = (int)(rx / dx); // Force to be integer
// Inject wavelet
for (int i = 0; i < ny; i++) {
// printf("i: %d\n", i);
UUo[irx][i] += dt * dt * slice[i][it];
}
// Define Courant numbers (squared)
double dtdx2 = (dt / dx) * (dt / dx);
double dtdy2 = (dt / dy) * (dt / dy);
// Update solution
for (int i = 4; i < nx - 4; i++) {
for (int j = 4; j < ny - 4; j++) {
UUm[i][j] = 2 * UUo[i][j] - UUm[i][j] + dtdx2 * v[i][j] * v[i][j] * (
-1.0/560 * UUo[i - 4][j]
+ 8.0/315 * UUo[i - 3][j]
- 1.0/5 * UUo[i - 2][j]
+ 8.0/5 * UUo[i - 1][j]
- 205.0/72 * UUo[i][j]
+ 8.0/5 * UUo[i + 1][j]
- 1.0/5 * UUo[i + 2][j]
+ 8.0/315 * UUo[i + 3][j]
- 1.0/560 * UUo[i + 4][j]) + dtdy2 * v[i][j] * v[i][j] * (
-1.0/560 * UUo[i][j - 4]
+ 8.0/315 * UUo[i][j - 3]
- 1.0/5 * UUo[i][j - 2]
+ 8.0/5 * UUo[i][j - 1]
- 205.0/72 * UUo[i][j]
+ 8.0/5 * UUo[i][j + 1]
- 1.0/5 * UUo[i][j + 2]
+ 8.0/315 * UUo[i][j + 3]
- 1.0/560 * UUo[i][j + 4]);
}
}
}
void awe_2d_explicit_solver_heterogeneous_8th_order(
double** UUo, double** UUm, int nx, int ny, double dx, double dy,
double dt, double** v, double* F, int it, double sx, double sy
) {
// Define Courant numbers (squared)
double dtdx2 = (dt / dx) * (dt / dx);
double dtdy2 = (dt / dy) * (dt / dy);
// Source location
int isx = (int)(sx / dx);
int isy = (int)(sy / dy);
// Inject wavelet
UUo[isx][isy] += dt * dt * F[it];
// Update solution without using slicing
for (int i = 4; i < nx - 4; i++) {
for (int j = 4; j < ny - 4; j++) {
UUm[i][j] = 2 * UUo[i][j] - UUm[i][j] + dtdx2 * v[i][j] * v[i][j] * (
-1.0/560 * UUo[i - 4][j]
+ 8.0/315 * UUo[i - 3][j]
- 1.0/5 * UUo[i - 2][j]
+ 8.0/5 * UUo[i - 1][j]
- 205.0/72 * UUo[i][j]
+ 8.0/5 * UUo[i + 1][j]
- 1.0/5 * UUo[i + 2][j]
+ 8.0/315 * UUo[i + 3][j]
- 1.0/560 * UUo[i + 4][j]) + dtdy2 * v[i][j] * v[i][j] * (
-1.0/560 * UUo[i][j - 4]
+ 8.0/315 * UUo[i][j - 3]
- 1.0/5 * UUo[i][j - 2]
+ 8.0/5 * UUo[i][j - 1]
- 205.0/72 * UUo[i][j]
+ 8.0/5 * UUo[i][j + 1]
- 1.0/5 * UUo[i][j + 2]
+ 8.0/315 * UUo[i][j + 3]
- 1.0/560 * UUo[i][j + 4]);
}
}
}
void append_2d_slice(double*** dest, double** src, int nx, int ny, int nt) {
for (int i = 0; i < nx; i++) {
for (int j = 0; j < ny; j++) {
dest[i][j][nt] = src[i][j]; // Debugger Points here Segmenation fault
}
}
}
void write_text_file_slice(double** slice, int nx, int ny, const char* filename) {
FILE *file = fopen(filename, "w");
if (file == NULL) {
fprintf(stderr, "Error opening file for writing.\n");
exit(1);
}
// Write slice data to the text file
for (int i = 0; i < nx; i++) {
for (int j = 0; j < ny; j++) {
fprintf(file, "%.50lf ", slice[i][j]);
}
fprintf(file, "\n");
}
fclose(file);
}
int main() {
// Define constants and arrays
int nz = 300;
int nx = 298;
double dz = 5.0;
double dx = 5.0;
double t0 = 0.05;
double CC = 0.5;
int nt = 2000;
double F[nt];
double sx = 25;
double sy = 800;
// Initialize UUo and UUm with zeros
double **UUo = malloc(nz * sizeof *UUo);
double **UUm = malloc(nz * sizeof *UUm);
double **UUo1 = malloc(nz * sizeof *UUo1);
double **UUm1 = malloc(nz * sizeof *UUm1);
double **v = malloc(nz * sizeof *v);
for (int i = 0; i < nz; i++) {
UUo[i] = calloc(nx, sizeof *UUo[i]);
UUm[i] = calloc(nx, sizeof *UUm[i]);
UUo1[i] = calloc(nx, sizeof *UUo1[i]);//addition
UUm1[i] = calloc(nx, sizeof *UUm1[i]);//addition
v[i] = calloc(nx, sizeof *v[i]);
}
// Initialize fff with zeros
double ***fff = malloc(nz * sizeof *fff);
for (int i = 0; i < nz; i++) {
fff[i] = calloc(nx, sizeof *fff[i]);
for (int j = 0; j < nx; j++) {
fff[i][j] = calloc(nt, sizeof *fff[i][j]);
}
}
// Initialize mmm with zeros
double ***mmm = malloc(nz * sizeof *mmm);
for (int i = 0; i < nz; i++) {
mmm[i] = calloc(nx, sizeof *mmm[i]);
for (int j = 0; j < nx; j++) {
mmm[i][j] = calloc(nt, sizeof *mmm[i][j]);
}
}
// Time stepping parameters
double dt = CC * dz / 5560.328125; // Assuming 'v' is the maximum value
double t[nt];
for (int i = 0; i < nt; i++) {
t[i] = i * dt;
}
// Define forcing function
double ss = 0.01;
for (int i = 0; i < nt; i++) {
F[i] = (1.0 - ((t[i] - t0) / ss) * ((t[i] - t0) / ss)) * exp(-((t[i] - t0) * (t[i] - t0)) / (2 * ss * ss));
}
// Read 'v' from the file
FILE *file = fopen("marmousi.txt", "r");
if (file == NULL) {
fprintf(stderr, "Error opening file.\n");
perror("fopen");
return 1; // Exit with an error code
}
for (int i = 0; i < nz; i++) {
for (int j = 0; j < nx; j++) {
if (fscanf(file, "%lf", &v[i][j]) != 1) {
fprintf(stderr, "Error reading from file.\n");
fclose(file);
return 1; // Exit with an error code
}
}
}
fclose(file);
printf("Forward Propagation...\n");
for (int it = 0; it < nt; it++) {
awe_2d_explicit_solver_heterogeneous_8th_order(UUo1, UUm1, nz, nx, dz, dx, dt, v, F, it, sx, sy);
// Swap UUo and UUm
double** tmp = UUo1;
UUo1 = UUm1;
UUm1 = tmp;
append_2d_slice(fff, UUo1, nz, nx, it);
}
printf("Forward Propagation Done !\n");
double** slice = fff[5];
write_text_file_slice(slice, nx, nt, "slice.txt");
printf("Backward Propagation...\n");
for (int it = nt - 1; it > 0; it--) {
//printf("%d\n",it);
awe_2d_heterogeneous_8th_order_data_time_reverse(UUo, UUm, nz, nx, dz, dx, dt, v, slice, it, sx);
double** tmp = UUo;
UUo = UUm;
UUm = tmp;
append_2d_slice(mmm, UUo, nz, nx, nt - 1 - it);// When I uncomment this part the code stops working and Debugger Points at line 93 on Windows Machine
}
printf("Backward Propagation Done !\n");
double** slicep = mmm[5];
write_text_file_slice(slicep, nx, nt, "slicep.txt");
// Free allocated memory
for (int i = 0; i < nz; i++) {
free(UUo[i]);
free(UUm[i]);
free(UUo1[i]);
free(UUm1[i]);
free(v[i]);
}
free(UUo);
free(UUm);
free(UUo1);
free(UUm1);
free(v);
for (int i = 0; i < nz; i++) {
for (int j = 0; j < nx; j++) {
free(fff[i][j]);
free(mmm[i][j]);
}
free(fff[i]);
free(mmm[i]);
}
free(fff);
free(mmm);
return 0;
}
这是脚本中使用的txt文件的链接:https://drive.google.com/drive/folders/1pmSkKrxS7IWDGyPiyl0SuNEh3gidmyFS?usp=drive_link
编辑:
上面带有未注释行
append_2d_slice(mmm, UUo, nz, nx, nt - 1 - it);
的代码适用于 Colab 和 wsl(适用于 Linux 的 Windows 子系统),但它不适用于 Windows 。问题是为什么?难道是gcc的版本?
“怀疑变量到函数的错误传递...”在这堆代码中很难分辨。虽然复制/粘贴是冲向终点线时的权宜解决方案,就像“龟兔赛跑”寓言一样,但权宜之计最终会自行决定。
下面是新产生的处理核心的第二个版本的简化(未经测试)。虽然确实应该避免过早的优化,但必须通过认识到人们何时走向复制功能的绝望深渊来调整这一经验法则。
void run(
double** UUo, double** UUm, double** v,
int nx, int ny,
double dx, double dy, double dt
) {
// coefficients as array of constants.
const double arr[] = { 205.0/72.0, 8.0/5.0, 1.0/5.0, 8.0/315.0, 1.0/560.0 };
// Courant numbers (squared)
const double dtdx2 = (dt / dx) * (dt / dx);
const double dtdy2 = (dt / dy) * (dt / dy);
for (int i = 4; i < nx - 4; i++)
for (int j = 4; j < ny - 4; j++) {
double someSqr = v[i][j] * v[i][j];
UUm[i][j]
= 2 * UUo[i][j] - UUm[i][j]
+ dtdx2 * someSqr * (
- arr[4] * UUo[i - 4][j]
+ arr[3] * UUo[i - 3][j]
- arr[2] * UUo[i - 2][j]
+ arr[1] * UUo[i - 1][j]
- arr[0] * UUo[i + 0][j]
+ arr[1] * UUo[i + 1][j]
- arr[2] * UUo[i + 2][j]
+ arr[3] * UUo[i + 3][j]
- arr[4] * UUo[i + 4][j] )
+ dtdy2 * someSqr * (
- arr[4] * UUo[i][j - 4]
+ arr[3] * UUo[i][j - 3]
- arr[2] * UUo[i][j - 2]
+ arr[1] * UUo[i][j - 1]
- arr[0] * UUo[i][j + 0]
+ arr[1] * UUo[i][j + 1]
- arr[2] * UUo[i][j + 2]
+ arr[3] * UUo[i][j + 3]
- arr[4] * UUo[i][j + 4] )
;
}
}
// Preliminary processing for this version...
void awe_2d_heterogeneous_8th_order_data_time_reverse(
double** UUo, double** UUm, int nx, int ny, double dx, double dy,
double dt, double** v, double** slice, int it, double rx
) {
// Get dimensions of the wavefield
int irx = (int)(rx / dx); // Force to be integer
// Inject wavelet
for (int i = 0; i < ny; i++)
UUo[irx][i] += dt * dt * slice[i][it];
// Update solution
run( UU0, UUm, v, nx, ny, dx, dy, dt );
}
// ... and preliminary processing for this version.
void awe_2d_explicit_solver_heterogeneous_8th_order(
double** UUo, double** UUm, int nx, int ny, double dx, double dy,
double dt, double** v, double* F, int it, double sx, double sy
) {
// Source location
int isx = (int)(sx / dx);
int isy = (int)(sy / dy);
// Inject wavelet
UUo[isx][isy] += dt * dt * F[it];
// Update solution without using slicing
run( UU0, UUm, v, nx, ny, dx, dy, dt );
}
虽然这可能会也可能不会回答OP的问题,但这种简化至少打开了更简单地追踪问题的窗口。 (第二个版本中单个元素的设置确实看起来很奇怪。是否缺少循环?)