プログラミング学習の記録 #22（C）：拡散

ボックスモデルにおける拡散方程式において、各ボックスで拡散係数が異なる場合のプログラムコードを以下に記録しておこうと思う。プログラミング言語は、C言語を用いている。

拡散方程式

$$
\dfrac{\partial D}{\partial t}
= y \cdot \dfrac{\partial}{\partial z} \left( \dfrac{\partial D}{\partial z} \right)
$$

各ボックスで拡散係数が異なる場合の拡散

$$
\Delta D_{m}
= \dfrac{y_{m-1} + y_{m}}{2} \cdot \dfrac{D_{m-1} - D_{m}}{(z_{m} - z_{m-1})^{2}}
+ \dfrac{y_{m+1} + y_{m}}{2} \cdot \dfrac{D_{m+1} - D_{m}}{(z_{m+1} - z_{m})^{2}}
$$

拡散のプログラムコード

海洋の鉛直構造を想定している。第$${m}$$層の深さ$${z_{m}}$$として、各ボックスの拡散係数$${y_{m}(=yD_{m})}$$とし、物理量$${D}$$の鉛直方向の拡散を考える。

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// Input the numbers of layer and depth.
#define layer 100
#define depth 600.0
// Input the nunmbers of time steps and end time.
#define year_rate 0.0027397300/2.0000
#define year_max 10.0

// the function of diffusion.
double get_diffusion(int m,double y0,double y,double y1,double M0,double M,double M1,double z0,double z,double z1)
{
    double diff;
    diff = ((y0 + y) / 2.0) * (M0 - M) / pow((z - z0),2) + ((y1 + y) / 2.0) * (M1 - M) / pow((z1 - z),2) ;
    return diff;
}

int main()
{
    FILE *fp;
    int m;
    double t,dt;
    double D[layer],diff_D[layer];
    double y[layer],yD[layer];
    double z[layer],z_bottom[layer],dz[layer],z_surface;

    // reseting time.
    t = 0.0 ;
    dt = year_rate ;

    // depth.
    z_surface = 0.0 ;
    z_bottom[0] = depth / (1.0 * layer) ;
    dz[0] = z_bottom[0] - z_surface ;
    z[0] = (z_bottom[0] + z_surface) / 2.0 ;
    for(m = 1 ; m < layer ; m = m + 1)
    {
        z_bottom[m] = z_bottom[m-1] + depth / (1.0 * layer) ;
        dz[m] = z_bottom[m] - z_bottom[m-1] ;
        z[m] = (z_bottom[m] + z_bottom[m-1]) / 2.0 ;
    }

    // initial numbers of D.
    for(m = 0 ; m < layer ; m = m + 1)
    {
        D[m] = 0.0 ;
    }

    // Input initial conditions.
    D[1] = 400.0 ;
    // OK.

    // Input conditions of diffusion coefficient.
    for(m = 0 ; m < layer ; m = m + 1)
    {
        if(z[m] < 50.0) // surface mixed layer.
        {
            y[m] = 20000.0 ;
        }else if(z[m] < 200.0) // thermal layer.
        {
            y[m] = 5000.0 ;
        }else // deep layer.
        {
            y[m] = 2500.0 ;
        }
    }
    // OK.

    for(m = 0 ; m < layer ; m = m + 1)
    {
        yD[m] = y[m] ;
    }

    // Calculating time variation.
    while(t < year_max)
    {
        for(m = 0 ; m < layer ; m = m + 1)
        {
            if(m == 0)
            {
                diff_D[m] = get_diffusion(m,yD[m],yD[m],yD[m+1],D[m],D[m],D[m+1],2*z[m]-z[m+1],z[m],z[m+1]) ;
            }else if(m == layer - 1)
            {
                diff_D[m] = get_diffusion(m,yD[m-1],yD[m],yD[m],D[m-1],D[m],D[m],z[m-1],z[m],2*z[m]-z[m-1]) ;
            }else
            {
                diff_D[m] = get_diffusion(m,yD[m-1],yD[m],yD[m+1],D[m-1],D[m],D[m+1],z[m-1],z[m],z[m+1]) ;
            }
        }


        for(m = 0 ; m < layer ; m = m + 1)
        {
            D[m] = D[m] + diff_D[m] * dt ;
        }

        t = t + dt ;
    }

    // Outputing CSV.
    fp = fopen("test.csv","w") ;
    if(fp == NULL)
    {
        printf("File error.\n");
        exit(1);
    }
    fprintf(fp,"layer,depth,D\n");
    for(m = 0 ; m < layer ; m = m + 1)
    {
        fprintf(fp,"%2d,%5.1f,%10.5f\n",m,z[m],D[m]);
    }
    fclose(fp);

    return 0;
}

-----

動け！タイムライン