理研 並列化のサマースクール

• RIKEN AICS HPC Summer School 2012
  • 期間：平成24年 8月6日〜10日
  • 場所：理研　計算科学研究機構
  • 内容：Fortran・Cを使った並列化（MPI・OpenMP）

• • 自分のモチベーションとしては、シミュレーションなどを通じて”並列化”という言葉は知っていたが、全くの無勉強であり、講習会を導入として習おうということ

• まとめ（修了テスト）
• 用語集
  • 並列計算
    • OpenMP
    • MPI
  • LINPACK
  • ムーアの法則
  • アムダールの法則
  • グスタフソンの法則
  • TLB (translation lookaside buffer)
  • ベクトル計算
  • マルチコア
  • メモリバンド幅
  • スレッド
  • ノイマンズボトルネック
  • キャッシュ・キャッシュメモリ

平衡点を求める

• 状態の変化を表した式から平衡状態を求めたい
• 1次元のブリュッセレータで試してみる
• パッケージrootSolveの中のmultiroot関数を使う
  • multirootの中身でもRの逆行列（solve）をつかっているようだ
• 黒がX、赤がYの濃度
  • 解が安定かどうかはたぶんまた別のはなし
  • 解の一意性も別のはなし
  • 初期条件をいろいろ変えてみる必要がありそう
• 計算の時間の刻み幅（dt）の値で結果が変わってしまうので拡散係数の修正が必要
• こちらのような数値計算で方程式の解をだせるようにしたい

library(rootSolve)

# ブリュッセレータのパラメータ
A<-2
B<-4.17
X0<-c(A,B/A) # 平衡となるXとYの値

# シミュレーションのパラメータ
L<-1
dr<-0.01
Nr<-floor(L/dr)+1
dt<-0.00001

# 拡散係数
Dx<-0.0016
Dy<-0.008

# 各点のX,Yの濃度、初期条件
C<-matrix(X0,nr=2,nc=Nr)
Nint<-2
C<-C+c(sin(Nint*pi*1:Nr/(L*Nr)),sin(Nint*pi*1:Nr/(L*Nr))) # 初期条件

C[,1]<-C[,Nr]<-X0  #固定端

par(mfrow=c(1,2))
plot(seq(from=0,to=1,length=Nr),C[1,],type="l",xlab="r",ylab="X conc.",ylim=range(C),main="initial concentration")
par(new=T)
plot(seq(from=0,to=1,length=Nr),C[2,],type="l",xlab="r",ylab="Y conc.",ylim=range(C),main="initial concentration",col=2)


# ブリュッセレータ
f<-function(X){
	Y<-c(0,0)
	Y[1] <- A+X[1]^2*X[2]-B*X[1]-X[1]
	Y[2] <- B*X[1]-X[1]^2*X[2]
	return(Y)
	}

# ブリュッセレータと拡散
h<-function(C){
	
	C<-matrix(C,nr=2,nc=Nr)
	C0<-C
	# ブリュッセレータの計算
	for(l in 1:Nr)
	{
	X<-C[,l]
	dX<-f(X)*dt
	C[,l]<-X+dX
	}
	
	# 拡散の計算
	dC_r<-C[,-1]-C[,-Nr]
	Ir<- -c(Dx,Dy)*dC_r
	C[,-1]<-C[,-1]+Ir
	C[,-Nr]<-C[,-Nr]-Ir
	
	# 固定端
	C[,1]<-C[,Nr]<-X0
	return(C-C0)	# 変化量をかえす
	}

# パッケージrootSolveの中のmultirootを使う
root<-multiroot(h,C,maxiter=100)$root
# print(matrix(h(root),nr=2,nc=Nr))

CC<-matrix(root,nr=2,nc=Nr)
plot(seq(from=0,to=1,length=Nr),CC[1,],type="l",xlab="r",ylab="X conc.",ylim=range(CC),main=paste(" equilibirium"))
par(new=T)
plot(seq(from=0,to=1,length=Nr),CC[2,],type="l",xlab="r",ylab="Y conc.",ylim=range(CC),main=paste(" equilibirium"),col=2)
# plot(seq(from=0,to=1,length=Nr),h(CC)[1,],type="l",ylim=c(-0.1,0.1),xlab="r",ylab="X conc.",main=paste("flux"))

• 下はシミュレーション用のコード

A<-2
B<-4.17
X0<-c(A,B/A) # 平衡となるXとYの値

L<-1
dr<-0.01
Nr<-floor(L/dr)+1

T<-3
dt<-0.00001
Nt<-T/dt

# 拡散係数
Dx<-0.0016
Dy<-0.008

# 各点のX,Yの濃度
C<-matrix(X0,nrow=2,ncol=Nr)
Nint<-2
C<-C+c(sin(Nint*pi*1:Nr/(L*Nr)),sin(Nint*pi*1:Nr/(L*Nr)))

range<-range(C)
plot(seq(from=0,to=1,length=Nr),C[1,],type="l",ylim=range,xlab="r",ylab="X conc.",main="time = 0")
par(new=T)
plot(seq(from=0,to=1,length=Nr),C[2,],type="l",ylim=range,xlab="r",ylab="Y conc.",main="time = 0",col=2)

# ブリュッセレータ
f<-function(X){
	Y<-c(0,0)
	Y[1] <- A+X[1]^2*X[2]-B*X[1]-X[1]
	Y[2] <- B*X[1]-X[1]^2*X[2]
	return(Y)
	}


for(t in 1:Nt)
{	
	# ブリュッセレータの計算
	for(l in 1:Nr)
	{
	X<-C[,l]
	dX<-f(X)*dt
	C[,l]<-X+dX
	}
	
	# 拡散の計算
	dC_r<-C[,-1]-C[,-Nr]
	Ir<- -c(Dx,Dy)*dC_r
	C[,-1]<-C[,-1]+Ir
	C[,-Nr]<-C[,-Nr]-Ir
	
	# 固定端
	C[,1]<-C[,Nr]<-X0

	if(t %% 100 == 0){
	plot(seq(from=0,to=1,length=Nr),C[1,],type="l",ylim=range,xlab="r",ylab="X conc.",main=paste("time = ",t*dt,sep=""))
	par(new=T)
	plot(seq(from=0,to=1,length=Nr),C[2,],type="l",ylim=range,xlab="r",ylab="Y conc.",main=paste("time = ",t*dt,sep=""),col=2)
	}

}

• こちらとこちらのpdf

library(fields)

# 連立方程式
f<-function(x){
z<-rep(0,dim)
z[1]<-(x[1]-3)^2+x[2]^2-3
z[2]<-sin(x[1])+exp(x[2]-1)-1
return(z)
}

# 偏微分係数を近似計算するときに差分を与える関数
g<-function(i){
	v<-rep(0,dim)
	v[i]<-1
	return(v)
	}



d<-0.001 # 偏微分係数の近似のための差分
dim<-2
A<-matrix(0,nr=dim,nc=dim)  #ヤコビアン

vectors<-c()
errors<-c()

# 初期値の集合 
by<-0.05
x0<-seq(from=0,to=5,by=by)
y0<-seq(from=-2.5,to=2.5,by=by)

points<-data.matrix(expand.grid(x0,y0))

lx<-length(x0)
ly<-length(y0)
L<-nrow(points)

pmat_x<-matrix(points[,1],lx,ly)
pmat_y<-matrix(points[,2],lx,ly)

vectors<-points

for(l in 1:L){
x<-XY<-points[l,]

# ヤコビアンの計算
for(i in 1:dim){
for(j in 1:dim){
	A[i,j]<-(f(x+d*g(j))[i]-f(x-d*g(j))[i])/(2*d)
}}

dXY<- try(-solve(A)%*%f(x))

if(is(dXY)[1] != "try-error"){# 連立方程式が解けた時
	vectors[l,]<-dXY
	}else{  # 連立方程式が解けなかった時
	errors<-cbind(errors,XY)
	vectors[l,]<-c(0,0)  # 0 を入れておく
	}

}

# ベクトルでプロット
#png("vector field.png")
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),type="n")

arrow.plot(points[,1],points[,2],vectors[,1],vectors[,2],arrow.ex=.2, length=.1,col=tim.colors(5)[(atan2(vectors[,2],vectors[,1])+pi)/(2*pi)*4+1], lwd=2,xlim=range(points[1,]),ylim=range(points[2,])) 
# dev.off()


# 境界を求める
vmat_x<-matrix(vectors[,1],lx,ly)
vmat_y<-matrix(vectors[,2],lx,ly)

bound_x<-(vmat_x[-1,]*vmat_x[-lx,]<0)
bound_y<-(vmat_y[,-1]*vmat_y[,-ly]<0)

# xの境界
bx_x<-(pmat_x[-1,][bound_x == 1]+pmat_x[-lx,][bound_x == 1])/2
bx_y<-(pmat_y[-1,][bound_x == 1]+pmat_y[-lx,][bound_x == 1])/2

# yの境界
by_x<-(pmat_x[,-1][bound_y == 1]+pmat_x[,-ly][bound_y == 1])/2
by_y<-(pmat_y[,-1][bound_y == 1]+pmat_y[,-ly][bound_y == 1])/2



# どっち？
# vmat or vmat * pmat ?? 

xstable<-(vmat_x[-1,]<0)*(vmat_x[-lx,]>0)
xstable_x<-(pmat_x[-1,][xstable == 1]+pmat_x[-lx,][xstable == 1])/2
xstable_y<-(pmat_y[-1,][xstable == 1]+pmat_y[-lx,][xstable == 1])/2

xunstable<-(vmat_x[-1,]>0)*(vmat_x[-lx,]<0)
xunstable_x<-(pmat_x[-1,][xunstable == 1]+pmat_x[-lx,][xunstable == 1])/2
xunstable_y<-(pmat_y[-1,][xunstable == 1]+pmat_y[-lx,][xunstable == 1])/2

ystable<-(vmat_y[,-1]<0)*(vmat_y[,-ly]>0)
ystable_x<-(pmat_x[,-1][ystable == 1]+pmat_x[,-ly][ystable == 1])/2
ystable_y<-(pmat_y[,-1][ystable == 1]+pmat_y[,-ly][ystable == 1])/2

yunstable<-(vmat_y[,-1]>0)*(vmat_y[,-ly]<0)
yunstable_x<-(pmat_x[,-1][yunstable == 1]+pmat_x[,-ly][yunstable == 1])/2
yunstable_y<-(pmat_y[,-1][yunstable == 1]+pmat_y[,-ly][yunstable == 1])/2








# 角度と境界をプロット
dev.new()
# png("angle.png")
par(mfrow=c(2,2))

#1
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),col=tim.colors(501)[(atan2(vectors[,2],vectors[,1])+pi)/(2*pi)*500+1], lwd=2,pch=15,cex=1,main="angle",xlab="X",ylab="Y") 

#2
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),col=tim.colors(501)[(atan2(vectors[,2],vectors[,1])+pi)/(2*pi)*500+1], lwd=2,pch=15,cex=1,main="merge",xlab="X",ylab="Y") 

par(new=T)
plot(xstable_x,xstable_y,pch=16,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(xunstable_x,xunstable_y,pch=3,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(ystable_x,ystable_y,pch=16,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(yunstable_x,yunstable_y,pch=3,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")

#3
plot(xstable_x,xstable_y,pch=16,xlim=range(points[,1]),ylim=range(points[,2]),main="boundary of x",xlab="X",ylab="Y")
par(new=T)
plot(xunstable_x,xunstable_y,pch=3,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="",)


#4
plot(ystable_x,ystable_y,pch=16,col="red",xlim=range(points[,1]),ylim=range(points[,2]),main="boundary of y",xlab="X",ylab="Y")
par(new=T)
plot(yunstable_x,yunstable_y,pch=3,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")

# dev.off()

# 安定性
dev.new()
# png("Stability.png",width=960)
par(mfrow=c(1,2))

#1
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),col=tim.colors(501)[ifelse((vmat_x<0),100,500)], lwd=2,pch=15,cex=1,main="x stability",xlab="X",ylab="Y")

par(new=T)
plot(xstable_x,xstable_y,pch=16,xlim=range(points[,1]),ylim=range(points[,2]),main="",xlab="",ylab="")
par(new=T)
plot(xunstable_x,xunstable_y,pch=3,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")

#2
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),col=tim.colors(501)[ifelse((vmat_y<0),100,500)], lwd=2,pch=15,cex=1,main="y stability",xlab="X",ylab="Y")

par(new=T)
plot(ystable_x,ystable_y,pch=16,col="red",xlim=range(points[,1]),ylim=range(points[,2]),main="",xlab="",ylab="")
par(new=T)
plot(yunstable_x,yunstable_y,pch=3,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
# dev.off()

• こちらに追加

数値計算

• 連立方程式の計算
• 一般化、大規模化のためにこちらの改良
  • Rの逆行列を用いている

x<-c(0,0)
T<-100
d<-0.01

dim<-length(x)

f<-function(x){
z<-rep(0,dim)
z[1]<-(x[1]-3)^2+x[2]^2-3
z[2]<-sin(x[1])+exp(x[2]-1)-1
return(z)
}


g<-function(i){
	v<-rep(0,dim)
	v[i]<-1
	return(v)
	}

A<-matrix(0,nr=dim,nc=dim)


for(t in 1:T){

for(i in 1:dim){
for(j in 1:dim){
	A[i,j]<-(f(x+d*g(j))[i]-f(x-d*g(j))[i])/(2*d)
}}


x<-x-solve(A)%*%f(x)

}

print(x)
print(f(x))

RとCを使う改良版

• RとCを使うときのファイル操作を、こちらのようにバイナリで行う
• Cのスクリプト
  • バイナリファイルを書きだす
  • 変数Aを書きだす部分をループの外に出した

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

#define PI 3.1415926535

int main (int argc, char *argv[]){
	
	int Nx, Ny, Nz;
	int i,j,k;
	
	FILE *fp;
	
	Nx = atoi(argv[1]);
	Ny = atoi(argv[2]);
	Nz = atoi(argv[3]);
	
//	printf("argv = %d\n",atoi(argv[1]));
	
	double A[Nx][Ny][Nz];
	
	for (k = 0; k < Nz; k++) {
	for (j = 0; j < Ny; j++) {
	for (i = 0; i < Nx; i++) {
		
		A[i][j][k] = (i + 1)*(j + 1)*(k + 1)*1.0;
		
//		fp = fopen("Numbers.txt","a");
//		if(fp == NULL)return 1;
//		fprintf(fp,"%f\n",A[i][j][k]);
//		fclose(fp);
		
	}}}
	
	
	fp = fopen("Numbers.txt","wb");
	if(fp == NULL){
		printf("ファイル操作中にエラー");
		exit(1);
	}
	fwrite(A,sizeof(double),Nx*Ny*Nz,fp);
	fclose(fp);
	
	return 0;
}

• Rのスクリプト

# mat.R

library(rgl)

Nx <- as.numeric(commandArgs(trailingOnly = TRUE)[1])
Ny <- as.numeric(commandArgs(trailingOnly = TRUE)[2])
Nz <- as.numeric(commandArgs(trailingOnly = TRUE)[3])

# m <- data.matrix(read.table("Numbers.txt",header=FALSE,sep=""))
m <-readBin("Numbers.txt",double(),n=Nx*Ny*Nz,endian="little")

A<-array(m,c(Nx,Ny,Nz))

plot3d(slice.index(A,1),slice.index(A,2),slice.index(A,3),col=rainbow(256)[A/(max(A))*256+1],alpha=ifelse(A!=0,1,0))
rgl.viewpoint(theta=30,phi=30)
rgl.snapshot("Array.png")

• シェルスクリプトはほぼそのまま

#!/bin/sh

Nx=30
Ny=30
Nz=30

gcc -o mat mat.c
./mat $Nx $Ny $Nz

R --vanilla --slave --args < mat.R $Nx $Ny $Nz

open Array.png
# rm Numbers.txt

バイナリで読み書き

• こちらでRの中でC++の関数を使えるようにしている
• そのやり方はわからないが、このときにCとRを使うことを考えていたのでもう少し改良
• CとRの間のファイルのやり取りをバイナリで行う
  • Rのバイナリファイル
  • Cについてはこちらのなかのバイナリデータ
    • Cについての以前の記事

• Cでバイナリファイルの読み書き
  • double型の変数が6個
  • バイナリで書かれたtest.txtも出力される

// mat.c

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

int main(void)
{
	double x[]={0.0,1.0,2.0,3.0,4.0,5.0} , y[6];

	FILE *fp;
	fp = fopen("test.txt","w");
	if(fp == NULL){
		printf("ファイル操作中にエラー");
		exit(1);
	}
	fwrite(x,sizeof(double),6,fp);
	fclose(fp);
		   
	fp = fopen("test.txt","r");
	if(fp == NULL){
		printf("ファイル操作中にエラー");
		exit(1);
	}
	fread(y,sizeof(double),6,fp);
	fclose(fp);
	
// 和を出力させている
	printf("%lf",y[0]+y[1]+y[2]+y[3]+y[4]+y[5]);
	
	return 0;
}

• Rで読み込み
  • double型変数6個の読み込み
  • readBin()を使う

# ファイルの場所とバイナリモードの読み込みの指定
to.read = file("~/Desktop/read/test.txt","rb")

# 読み込み方の指定をする
readBin(to.read,double(),n=6,endian = "little")

• Rの実行結果

> to.read = file("~/Desktop/read/test.txt","rb")
> readBin(to.read,double(),n=6,endian = "little")
[1] 0 1 2 3 4 5

CとRを動かす

• 久しぶりにCの復習をしてみた
• シェルでCとRを動かせるようにシェルスクリプトを書いた
  • Cを計算に使いたい
  • Rで図を描きたい
• いままでマスで仕切ってマスごとに計算するというシミュレーションをやってきていたので、Cでも同じことをしてみる
• マスの大きさをシェルスクリプトから与えられるようにしてみた
  • Cのコマンドライン引数の使い方
  • Rのコマンドライン引数の使い方
• 下の3つを同じディレクトリに置いておいて、シェルスクリプトを実行すればよい
• 全体の流れとしては
  • シェルで3次元空間の大きさを与える（Nx、Ny、Nz）
  • CがNx、Ny、Nzの大きさの空間内の各マスの中身を計算する（3次元の配列）
  • Cが3次元の配列の中身を書き出す（Numbers.txt）
  • RがNumbers.txtを読み込む
  • Rのrglパッケージで3dplotして図を描きだす（Array.png）
  • シェルでArray.pngを開き、Numbers.txtは消去する

#!/bin/sh

Nx=30
Ny=30
Nz=30

# Cをコンパイル、実行まで
gcc -o mat mat.c
./mat $Nx $Ny $Nz

# Rを実行
R --vanilla --slave --args < mat.R $Nx $Ny $Nz

# ファイルの実行と消去
open Array.png
rm Numbers.txt

/* mat.c */

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

#define PI 3.1415926535

int main (int argc, char *argv[]){
	
	int Nx, Ny, Nz;
	int i,j,k;
	
	FILE *fpA;
	
	Nx = atoi(argv[1]);
	Ny = atoi(argv[2]);
	Nz = atoi(argv[3]);
		
	double A[Nx][Ny][Nz];
	
	for (k = 0; k < Nz; k++) {
	for (j = 0; j < Ny; j++) {
	for (i = 0; i < Nx; i++) {
		
		A[i][j][k] = (i + 1)*(j + 1)*(k + 1)*1.0;
		
		fpA = fopen("Numbers.txt","a");
		if(fpA == NULL)return 1;
		fprintf(fpA,"%f\n",A[i][j][k]);
		fclose(fpA);
	}}}
	
	return 0;
}

# mat.R
library(rgl)

Nx <- as.numeric(commandArgs(trailingOnly = TRUE)[1])
Ny <- as.numeric(commandArgs(trailingOnly = TRUE)[2])
Nz <- as.numeric(commandArgs(trailingOnly = TRUE)[3])

m <- data.matrix(read.table("Numbers.txt",header=FALSE,sep=""))
A<-array(m,c(Nx,Ny,Nz))

plot3d(slice.index(A,1),slice.index(A,2),slice.index(A,3),col=rainbow(256)[A/(max(A))*256+1],alpha=ifelse(A!=0,1,0))
rgl.snapshot("Array.png")

ローレンツアトラクタ

• ローレンツ方程式（Wiki）
• 位相次元
• フラクタル次元
  • フラクタル次元の中に、情報次元、ボックスカウント次元、ハウスドルフ次元、レニー次元などの定義があるらしい

# ローレンツ方程式
p<-10
r<-28
b<-8/3
f<-function(x,y,z){-p*x+p*y}
g<-function(x,y,z){-x*z+r*x-y}
h<-function(x,y,z){x*y-b*z}

# 平衡点を求める
# 初期の値
x<-10
y<-10
z<-10

T<-100
d<-0.01

for(t in 1:T){

A<-matrix(c(
(f(x+d,y,z)-f(x-d,y,z))/(2*d),(f(x,y+d,z)-f(x,y-d,z))/(2*d),(f(x,y,z+d)-f(x,y,z-d))/(2*d),
(g(x+d,y,z)-g(x-d,y,z))/(2*d),(g(x,y+d,z)-g(x,y-d,z))/(2*d),(g(x,y,z+d)-g(x,y,z-d))/(2*d),
(h(x+d,y,z)-h(x-d,y,z))/(2*d),(h(x,y+d,z)-h(x,y-d,z))/(2*d),(h(x,y,z+d)-h(x,y,z-d))/(2*d)
),byrow=T,3,3)

dX<- -solve(A)%*%c(f(x,y,z),g(x,y,z),h(x,y,z))
x<-x+dX[1]
y<-y+dX[2]
z<-z+dX[3]
}

x;y;z
f(x,y,z);g(x,y,z);h(x,y,z)  # 確認

# 軌跡を書いてみる
library(rgl)

x<-0.1;y<-0;z<-0
V<-matrix(c(x,y,z),3,1)
dt<-0.01
T<-100
Nt<-T/dt

for(t in 1:Nt){
dX<-c(f(x,y,z),g(x,y,z),h(x,y,z))*dt
x<-x+dX[1]
y<-y+dX[2]
z<-z+dX[3]
V<-cbind(V,c(x,y,z))
}

plot3d(t(V))

library(rgl)

# 平行移動と拡大縮小で(0,1)におさめる
U<-V
U[1,]<-(V[1,]-min(V[1,]))/(max(V[1,])-min(V[1,]))
U[2,]<-(V[2,]-min(V[2,]))/(max(V[2,])-min(V[2,]))
U[3,]<-(V[3,]-min(V[3,]))/(max(V[3,])-min(V[3,]))
plot3d(t(U),xlim=c(0,1),ylim=c(0,1),zlim=c(0,1))

# 精度を与えてデータを集計する
L<-40
M<-array(0,c(L,L,L))
for(l in 1:dim(U)[2]){
i<-floor(U[1,l]*L)
j<-floor(U[2,l]*L)
k<-floor(U[3,l]*L)
M[i,j,k]<-M[i,j,k]+1
}

plot3d(slice.index(M,1),slice.index(M,2),slice.index(M,3),alpha=ifelse(M==0,0,1))

# このままでは重いので点のあるところのみを選んでプロット
M2<-rbind((which(M > 0)-1) %% L,(which(M > 0)-1) %% L^2,(which(M > 0)-1) %% L^3)+1
plot3d(t(M2))

# 各位置での速度を計算する
Nx<-Ny<-Nz<-40

v1<-seq(from=min(V[1,],to=max(V[1,])),length=Nx)
v2<-seq(from=min(V[2,],to=max(V[2,])),length=Ny)
v3<-seq(from=min(V[3,],to=max(V[3,])),length=Nz)

# 各位置での速度を計算する
F<-array(0,c(Nx,Ny,Nz,3))
for(x in 1:Nx){
for(y in 1:Ny){
for(z in 1:Nz){
lx<-v1[x]
ly<-v2[y]
lz<-v3[z]
F[x,y,z,1]<-f(lx,ly,lz)
F[x,y,z,2]<-g(lx,ly,lz)
F[x,y,z,3]<-h(lx,ly,lz)
}}}

# 速度を色にしてプロットする
Col1<-(F[,,,1]-min(F[,,,1]))/(max(F[,,,1])-min(F[,,,1]))
Col2<-(F[,,,2]-min(F[,,,2]))/(max(F[,,,2])-min(F[,,,2]))
Col3<-(F[,,,3]-min(F[,,,3]))/(max(F[,,,3])-min(F[,,,3]))


plot3d(v1[slice.index(F,1)],v2[slice.index(F,2)],v3[slice.index(F,3)],xlab="x",ylab="y",zlab="z",col=rainbow(101)[Col1*100+1],main="velocity x")
plot3d(v1[slice.index(F,1)],v2[slice.index(F,2)],v3[slice.index(F,3)],xlab="x",ylab="y",zlab="z",col=rainbow(101)[Col2*100+1],main="velocity y")
plot3d(v1[slice.index(F,1)],v2[slice.index(F,2)],v3[slice.index(F,3)],xlab="x",ylab="y",zlab="z",col=rainbow(101)[Col3*100+1],main="velocity z")

# 偏微分係数を計算してプロット
# z方向について
Fz<-F[,,-1,3]-F[,,-Nz,3]

err<-0.000001
if(max(Fz) - min(Fz) < err){
ColFz<-1
}else{
ColFz<-(Fz-min(Fz))/(max(Fz)-min(Fz))
}

plot3d(v1[slice.index(Fz,1)],v2[slice.index(Fz,2)],v3[slice.index(Fz,3)],xlab="x",ylab="y",zlab="z",col=rainbow(101)[ColFz*100+1],main="Fz")

# 時系列データから速度を計算し色にしてプロット
du1<-U[1,-1]-U[1,-length(U[1,])]
du1<-(du1-min(du1))/(max(du1)-min(du1))

du2<-U[2,-1]-U[2,-length(U[2,])]
du2<-(du2-min(du2))/(max(du2)-min(du2))

du3<-U[3,-3]-U[3,-length(U[3,])]
du3<-(du3-min(du3))/(max(du3)-min(du3))

C1<-c(du1[1],du1*100+1)
C2<-c(du2[1],du2*100+1)
C3<-c(du3[1],du3*100+1)

plot3d(U[1,],U[2,],U[3,],col=rainbow(101)[C1])
plot3d(U[1,],U[2,],U[3,],col=rainbow(101)[C2])
plot3d(U[1,],U[2,],U[3,],col=rainbow(101)[C3])