熱伝導シミュレーション

下記のソースコードは，有限差分法によって熱伝導方程式をシミュレーションするためのプログラムです．

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>

struct simulator {
    int n;
    float *T1;
    float *T2;
    int tmax;
    float alpha;
    float w1;
    float w2;
    float dx;
    float dt;
};

static struct timespec diff(struct timespec start, struct timespec end)
{
    struct timespec ts;

    if (end.tv_nsec - start.tv_nsec < 0) {
        ts.tv_sec = end.tv_sec - start.tv_sec - 1;
        ts.tv_nsec = 1000000000 + end.tv_nsec - start.tv_nsec;
    } else {
        ts.tv_sec = end.tv_sec - start.tv_sec;
        ts.tv_nsec = end.tv_nsec - start.tv_nsec;
    }
    return ts;
}

static void write_result(struct simulator *sim, int step)
{
    char fname[32];
    sprintf(fname, "result_%d.csv", step);

    FILE *fd = fopen(fname, "w");

    fprintf(fd, "x,y,T\n");

    for (int y = 1; y <= sim->n; y++) {
        for (int x = 1; x <= sim->n; x++) {
            fprintf(fd, "%f,%f,%f\n", (x - 1) * sim->dx, (y - 1) * sim->dx,
                    sim->T1[x + y * (sim->n + 2)]);
        }
    }

    fclose(fd);
}

static void step(struct simulator *sim)
{
    for (int y = 1; y <= sim->n; y++) {
        for (int x = 1; x <= sim->n; x++) {
            sim->T2[x + y * (sim->n + 2)] =
                sim->w1 * sim->T1[x + y * (sim->n + 2)] +
                sim->w2 * (sim->T1[(x - 1) + y * (sim->n + 2)] +
                           sim->T1[x + (y - 1) * (sim->n + 2)] +
                           sim->T1[(x + 1) + y * (sim->n + 2)] +
                           sim->T1[x + (y + 1) * (sim->n + 2)]);
        }
    }

    float *tmp = sim->T2;
    sim->T2 = sim->T1;
    sim->T1 = tmp;
}

static void compute(struct simulator *sim)
{
    for (int i = 0; i < sim->tmax; i++) {
        step(sim);

        if (i % 10 == 0) {
            printf("Writing out step %d\n", i);
            write_result(sim, i);
        }
    }
}

static void print_settings(struct simulator *sim)
{
    printf("grid size:  %dx%d\n", sim->n, sim->n);
    printf("tmax:       %d\n", sim->tmax);
    printf("alpha:      %f\n", sim->alpha);
    printf("delta x:    %f\n", sim->dx);
    printf("delta t:    %f\n", sim->dt);
}

static void initialize(struct simulator *sim)
{
    for (int y = 0; y <= sim->n + 1; y++) {
        for (int x = 0; x <= sim->n + 1; x++) {
            float temp = 50.0f;

            if (x == 0 || y == 0) {
                temp = 0.0f;
            } else if (x == sim->n + 1 || y == sim->n + 1) {
                temp = 100.0f;
            }

            sim->T1[x + y * (sim->n + 2)] = temp;
            sim->T2[x + y * (sim->n + 2)] = temp;
        }
    }
}

int main(int argc, char *argv[])
{
    int ch;
    int n = 100, tmax = 100;

    while ((ch = getopt(argc, argv, "t:n:")) != -1) {
        switch (ch) {
        case 'n':
            n = atoi(optarg);
            break;
        case 't':
            tmax = atoi(optarg);
            break;
        case 'h':
        default:
            fprintf(stderr, "Usage: %s -n <size>\n", argv[0]);
            break;
        }
    }

    struct simulator sim;

    sim.n = n;
    sim.T1 = (float *)calloc((n + 2) * (n + 2), sizeof(float));
    sim.T2 = (float *)calloc((n + 2) * (n + 2), sizeof(float));
    sim.tmax = tmax;

    sim.alpha = 0.1f;
    sim.dx = 0.01f;
    sim.dt = 0.0001f;
    sim.w1 = 1.0f - (4.0f * sim.alpha * sim.dt) / (sim.dx * sim.dx);
    sim.w2 = (sim.alpha * sim.dt) / (sim.dx * sim.dx);

    print_settings(&sim);

    initialize(&sim);

    struct timespec start_ts, end_ts;
    clock_gettime(CLOCK_MONOTONIC, &start_ts);

    compute(&sim);

    clock_gettime(CLOCK_MONOTONIC, &end_ts);

    struct timespec elapsed_ts = diff(start_ts, end_ts);

    printf("elapsed:    %ld.%09ld [s]\n", elapsed_ts.tv_sec,
           elapsed_ts.tv_nsec);

    return 0;
}

やってみよう

• このプログラムをコンパイルし，実行してみましょう．
• ParaViewを用いて，出力されたCSVファイルを可視化してみましょう．
• 初期条件を変更し，シミュレーション結果がどのように変わるか調べてみましょう．