Cで周期境界を再現したい

前提・実現したいこと

砂紋のパターン形成のシミュレーションを行いたい。

発生している問題・エラーメッセージ

上記のシミュレーションを行うプログラムを作成しています。
プログラム上で周期境界を再現したいのですが、やり方がどうしても分かりません。
簡単に言うと、x=50,y=50の幅の平面で2次元配列h[50][50]を扱う際に
[-1]など中身が負の数になってしまうとき、[-1]→[49]として[-1]をもう一方の端として扱うようにしたいです。
C言語に慣れておらず拙い説明で申し訳ありませんが、ご回答いただけると嬉しいです。

該当のソースコード

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


int main(void)
{
	int x_scale, y_scale, step,ir,jr;
	x_scale = 50; //x方向長さ
	y_scale = 50; //y方向長さ
    step = 1; //計算回数

	double yuragi, q, b, D, side, slant,v,w;
	yuragi = 0.3; //
	q = 0.1; //跳躍により移動する砂の量
	b = 1.0; //跳躍距離のパラメータ
	D = 0.1; //転がりにより移動する砂の量のパラメータ(転がり出る砂の量=Dh)
	side = 2.0; //転がり計算の時に使うパラメータ
	slant = 1.0; //同上

    srand(time(NULL));

	int i,j,s, t, L, k, L_0;//L=L_0 + bh (跳躍距離)
	L_0 = 3; //風速に対応するパラメータ

	double h[50][50];  //[x,y]における砂の高さ(積載量)

	
	FILE*fp;
	fp=fopen("output.dat", "w");
	for (j = 0;j < y_scale;j++)
	{
		for (i = 0;i < x_scale;i++)
		{
			h[i][j] = 2.85 + (yuragi * rand() / (double)RAND_MAX); //各座標のhの初期値を2.85~3.15の間で乱数で取る
		}
	}


	for (t = 1; t <= step; t++) 
	{

		for (k = 1; k < 2500; k++) //跳躍計算
		{   
		v = (int)(50 * rand() / (double)RAND_MAX);
		w = (int)(50 * rand() / (double)RAND_MAX);
        ir = v;
        jr = w;                                           //

			if (h[ir][jr] > 0)
			{
				L = L_0 + (b * (int)h[(ir)][(jr)]);
				h[(ir)][(jr)] -= q;
				h[(ir + L)][(jr)] += q;
			}
			else
			{
				continue;
			}

		}
		for (j = 1;j < y_scale - 1;j++) {
			for (i = 1;i < x_scale - 1;i++) {
				if (h[i][j] > h[(i - 1)][(j)] &&
					h[i][j] > h[(i + 1)][(j)] && h[i][j] > h[(i)][(j - 1)] &&
					h[i][j] > h[i][(j + 1)] && h[i][j] > h[(i - 1)][(j - 1)] &&
					h[i][j] > h[i - 1][(j + 1)] && h[i][j] > h[(i + 1)][(j - 1)] &&
					h[i][j] > h[i + 1][(j + 1)])
				{
					h[(i - 1)][(j)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant));
					h[(i + 1)][(j)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant));
					h[(i)][(j - 1)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant));
					h[(i)][(j + 1)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant));
					h[(i - 1)][(j - 1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant));
					h[(i - 1)][(j + 1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant));
					h[(i + 1)][(j - 1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant));
					h[(i + 1)][(j + 1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant));
					h[i][j] -= D * h[i][j];
				}
				else {
					continue;
				}
			}

		}

		for (j = 0;j < y_scale;j++) {
            for(i = 0; i < x_scale; i++ ){
			printf("%f\n",h[i][j]);
			fprintf(fp, "%f\n",h[i][j]);
            }
		}
		
	}
	
	fclose(fp);
}

変更後のソースコード

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


int main(void)
{
	int x_scale, y_scale, step;
	x_scale = 50; //x方向長さ
	y_scale = 50; //y方向長さ
    step = 1; //計算回数

	double yuragi, q, b, D, side, slant,v,w;
	yuragi = 0.3; //
	q = 0.1; //跳躍により移動する砂の量
	b = 1.0; //跳躍距離のパラメータ
	D = 0.1; //転がりにより移動する砂の量のパラメータ(転がり出る砂の量=Dh)
	side = 2.0; //転がり計算の時に使うパラメータ
	slant = 1.0; //同上

    srand(time(NULL));

	int x,y,i,j,i_1,i_2,j_1,j_2,i_0,s, t, L, k, L_0;//L=L_0 + bh (跳躍距離)
	L_0 = 3; //風速に対応するパラメータ
    i = (x_scale + (x % x_scale))%(x_scale);
    j = (y_scale + (y % y_scale))%(y_scale);
    i_1 = (x_scale + ((x - 1) % x_scale))%(x_scale);
    i_2 = (x_scale + ((x + 1) % x_scale))%(x_scale);
    j_1 = (y_scale + ((y - 1) % y_scale))%(y_scale);
    j_2 = (y_scale + ((y + 1) % y_scale))%(y_scale);
    i_0 = (x_scale + ((x + L) % x_scale))%(x_scale);


	double h[50][50];  //[x,y]における砂の高さ(積載量)

	
	FILE*fp;
	fp=fopen("output.dat", "w");
	for (y = 0;y < y_scale;y++)
	{
		for (x = 0;x < x_scale;x++)
		{
			h[x][y] = 2.85 + (yuragi * rand() / (double)RAND_MAX); //各座標のhの初期値を2.85~3.15の間で乱数で取る
		}
	}


	for (t = 1; t <= step; t++) 
	{

		for (k = 1; k < 2500; k++) //跳躍計算
		{   
		v = (int)(50 * rand() / (double)RAND_MAX);
		w = (int)(50 * rand() / (double)RAND_MAX);
        x = v;
        y = w;                                           
        L = L_0 + (b * (int)h[(i)][(j)]);
            
			if (h[i][j] > 0)
			{
				h[(i)][(j)] -= q;
				h[(i_0)][(j)] += q;
			}
			else
			{
				continue;
			}

		}

		for (y = 0;y < y_scale ;y++) 
        {
			for (x = 0;x < x_scale ;x++) 
            {
				if (h[i][j] > h[(i_1)][(j)] &&
					h[i][j] > h[(i_2)][(j)] && h[i][j] > h[(i)][(j_1)] &&
					h[i][j] > h[i][(j_2)] && h[i][j] > h[(i_1)][(j_1)] &&
					h[i][j] > h[i_1][(j_2)] && h[i][j] > h[(i_2)][(j_1)] &&
					h[i][j] > h[i_2][(j_2)])
				{
					h[(i_1)][(j)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant));
					h[(i_2)][(j)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant));
					h[(i)][(j_1)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant));
					h[(i)][(j_2)] += D * h[i][j] * (1.0 / 4.0) * (side / (side + slant));
					h[(i_1)][(j_1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant));
					h[(i_1)][(j_2)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant));
					h[(i_2)][(j_1)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant));
					h[(i_2)][(j_2)] += D * h[i][j] * (1.0 / 4.0) * (slant / (side + slant));
					h[i][j] -= D * h[i][j];
				}
				else {
					continue;
				}
			}

		}

		for (y = 0;y < y_scale;y++) {
            for(x = 0; x < x_scale; x++ ){
			printf("%f\n",h[i][j]);
			fprintf(fp, "%f\n",h[i][j]);
            }
		}
		
	}
	
	fclose(fp);
}

出力

-0.039956
-0.039956
-0.039956
-0.039956
↑これがずっと続きます

補足情報（FW/ツールのバージョンなど）

Visual Studio Codeで作成しています。

y_waiwai

2021/11/06 11:43

提示のコードはどういう動作をするんでしょうか

Apexer

2021/11/06 11:47

動作はするのですが、出力されたファイルに何も書かれていない、といった感じです。

BeatStar

2021/11/06 12:30

あと、コードは"<code>"または"<コード>"のボタンがあるはずなので、それを押して出てくるヤツの中に書きましょう。今のままでは読みづらいです。質問は編集できるので編集しましょう。

jimbe

2021/11/06 18:23

> 動作はするのですが、出力されたファイルに何も書かれていない計算うんぬんでは無く、ファイルが空になってしまうという問題なのでしょうか。

Apexer

2021/11/07 05:29 編集

>>BeatStarさんすいません、編集しました。

Apexer

2021/11/07 06:39 編集

majuponiさんの回答を参考に、プログラムを編集しました。数値は出力されるようになったのですが、出力された数値がすべて同じで、しかも負の値になってしまいます。(0以上を想定していた）おそらく乱数の使い方に誤りがあるとおもうのですが、よくわかりません。。

jimbe

2021/11/07 06:40

rand() による乱数の求め方を間違っていないでしょうか。 rand() のドキュメントや、 rand() を含む計算の結果を printf で確認しては如何でしょう。

Apexer

2021/11/07 07:05

rand()を含む計算の結果を確認したところ、やはりすべて同じ数値になりました。

jimbe

2021/11/07 08:13 編集

> すべて同じ数値・・・もう一声・・・といいますか、つまり乱数を利用する計算が間違っている、rand() の利用方法が間違っているとは言えないでしょうか? [Cライブラリ関数-rand（）](https://www.tutorialspoint.com/c_standard_library/c_function_rand.htm )

行動規範の内容に同意します

回答2件

C
1i=(w+n%w)%w;

この式が何を意味するか、どうしてこれでよいか考えてみましょう。

投稿2021/11/06 14:15

majiponi

総合スコア1722

Apexer

2021/11/07 06:41

剰余で表す考え方は盲点でした！nがどの整数値をとってもiは0~49に収まるということですよね。ありがとうございます！

Apexer

2021/11/07 07:33

数値は出力されるようになったのですが、出力された数値がすべて同じで、しかも負の値になってしまいます。(0以上を想定していた）乱数の使い方に誤りがあるとおもうのですが、どうすればいいかわかりません。。

行動規範の内容に同意します

周期境界とやらは存じませんが、今回の場合、そのアルゴリズムについての質問では無い為回答します。

...[-1]など中身が負の数になってしまうとき、[-1]→[49]として[-1]...

つまり、『負の数でアクセスする場合、強制的に最大値にする』的な意味ですよね?
それなら単純にそう実装すればいいだけでは。

そもそもプログラミングというのは『こう書けばいい』というものではありません。
言うなれば『現実世界のシミュレーション』です。

たとえば、この周期境界のやつを実際に手作業でやってみてください。
そして、今回の質問である、『負の数なら~』も同様にやってみてください。

そうすると、『要素数に相当する変数が負の数の場合』-> 『要素数の変数を最大値にする』となるはずです。

たとえば、

int pos = -1;
int desk[10][10];
int max = 10; // 各要素の最大値

// posが-1なら
if( pos == -1 ){
    // 強制的に最大値にする
    pos = max;
}

// 取得する
int n = desk[0][pos];

的な感じで。

ただし、負の数っていうのは -1　だけじゃなくて -2, -3, -4, ... とあります。
それをすべてやるのは非効率。何か規則性があるはずです。

数直線を思い浮かべてください。実際にノートに書いてもいいです。
とにかく、負の数を出来る限りチェックつけましょう。-1, -2, -3, -4, ... と。

で、そうすると『常に0より小さい』ことがわかりますね。そう、『0より小さいかどうか』で判断できますね。0は…含みませんね。

ということは、 pos＜0 の条件式を満たす場合が負の数ですね。

さあ、これでわかりましたね。後は組み込んで、デバッグするだけです。

投稿2021/11/06 11:52

BeatStar

総合スコア4962