Class | Alge_Solv |
In: |
alge_solv.f90
|
代数演算を用いて偏微分方程式を解くモジュール
Subroutine : | |||
x(:) : | real, intent(in)
| ||
y(:) : | real, intent(in)
| ||
rho(size(x),size(y)) : | real, intent(in)
| ||
eps : | real, intent(in)
| ||
boundary : | character(4), intent(in)
| ||
psi(size(x),size(y)) : | real, intent(inout)
| ||
bound_opt(size(x),size(y)) : | real, intent(in), optional
| ||
a(size(x),size(y)) : | real, intent(in), optional
| ||
b(size(x),size(y)) : | real, intent(in), optional
| ||
c(size(x),size(y)) : | real, intent(in), optional
| ||
d(size(x),size(y)) : | real, intent(in), optional
| ||
e(size(x),size(y)) : | real, intent(in), optional
| ||
undef : | real, intent(in), optional
| ||
inner_bound(size(x),size(y)) : | integer, intent(in), optional
| ||
init_flag : | logical, intent(in), optional
|
ガウス=ザイデル法によるポアソン方程式の求積 各オプション配列は, ポアソン系の各微分項の値. デフォルトはゼロで設定される. $$a\dfrac{partial ^2\psi}{partial x^2} +b\dfrac{partial ^2\psi}{partial x\partial y} +c\dfrac{partial ^2\psi}{partial y^2} +d\dfrac{partial psi}{partial x} +e\dfrac{partial psi}{partial y} =rho $$ の各係数に対応している.
subroutine Poisson_GauSei(x, y, rho, eps, boundary, psi, bound_opt, a, b, c, d, e, undef, inner_bound, init_flag ) ! ガウス=ザイデル法によるポアソン方程式の求積 ! 各オプション配列は, ポアソン系の各微分項の値. デフォルトはゼロで設定される. ! $$a\dfrac{\partial ^2\psi}{\partial x^2} +b\dfrac{\partial ^2\psi}{\partial x\partial y} +c\dfrac{\partial ^2\psi}{\partial y^2} +d\dfrac{\partial \psi}{\partial x} +e\dfrac{\partial \psi}{\partial y} =\rho $$ ! の各係数に対応している. implicit none real, intent(in) :: x(:) ! 領域の横座標 real, intent(in) :: y(:) ! 領域の縦座標 real, intent(in) :: rho(size(x),size(y)) ! ポアソン方程式の強制項 ! rho =0 でラプラス方程式も求積可能 real, intent(in) :: eps ! 収束条件 character(4), intent(in) :: boundary ! 境界条件 ! 4 文字で各辺の境界条件を与える. ! 1 文字目 : x 下端, 2 文字目 : y 左端, 3 文字目 : x 上端, ! 4 文字目 : y 右端 ! boundary は 1 : 固定端境界, 2 : 自由端境界, 3 : 周期境界 real, intent(in), optional :: bound_opt(size(x),size(y)) ! 境界での強制 ! ノイマン境界の場合 : フラックス値 real, intent(in), optional :: a(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: b(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: c(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: d(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: e(size(x),size(y)) ! 各微分項の係数 real, intent(inout) :: psi(size(x),size(y)) ! ポアソン方程式の解 real, intent(in), optional :: undef ! 未定義値 integer, intent(in), optional :: inner_bound(size(x),size(y)) ! 内部領域の境界. 値に応じてその格子点で境界値計算 ! 1 = 固定端境界, 10 = 境界の内側. ! 2 = y 方向自由端境界 (フラックスは上向き) ! -2 = y 方向自由端境界 (フラックスは下向き) ! 4 = x 方向自由端境界 (フラックスは右向き) ! -4 = x 方向自由端境界 (フラックスは左向き) ! 3 = 周期境界, -3 = 隅領域で両方とも周期境界 ! 8 = |_, ~| で両方とも自由境界条件 ! -8 = |~, _| で両方とも自由境界条件 ! この引数が与えられなければ全領域を計算する. ! 境界の内側格子点 (10) は反復計算を行わず, ! undef で設定された値もしくはゼロが入る. ! このときの境界値は bound_opt の値が用いられる. logical, intent(in), optional :: init_flag ! psi の値をゼロで初期化するか. ! .true. = 初期化する. .false. = 初期化しない. ! デフォルトでは初期化する. integer :: i, j integer :: nx ! x 方向の配列要素 integer :: ny ! y 方向の配列要素 integer :: signb ! 各係数を計算するかどうか integer, dimension(size(x),size(y)) :: ib real :: defun real :: tmp, err, err_max real :: tmp_b real :: bnd(size(x),size(y)) real :: dx(size(x)), dy(size(y)), dx2(size(x)), dy2(size(y)) real, dimension(size(x),size(y)) :: dxdy real, dimension(size(x),size(y)) :: at, bt, ct, dt, et real, dimension(size(x),size(y)) :: adp, adm, cep, cem, ac, divi character(4) :: bound logical, dimension(size(x),size(y)) :: inner_flag bound(1:4)=boundary(1:4) nx=size(x) ny=size(y) !-- 応答関数の初期化 if(present(init_flag))then if(init_flag.eqv..true.)then psi = 0.0 end if else psi = 0.0 end if !-- 内部境界の判別フラグの設定 if(present(inner_bound))then call set_bound( bound, ib, inner_flag, inner_bound ) else call set_bound( bound, ib, inner_flag ) end if !-- 領域・内部境界における境界値の設定 if(present(bound_opt))then call setval_bound( ib, bnd, psi, bound_opt ) else call setval_bound( ib, bnd, psi ) end if !-- 未定義値の設定 if(present(undef))then defun=undef else defun=0.0 end if !-- 係数の代入 !-- a, c については, 値が入れられていなければ, 全配列に 1 を代入する. if(present(a))then call set_coe( at, ext=a ) else call set_coe( at, def=1.0 ) end if if(present(c))then call set_coe( ct, ext=c ) else call set_coe( ct, def=1.0 ) end if if(present(b))then call set_coe( bt, ext=b ) signb=1 else call set_coe( bt, def=0.0 ) signb=0 end if if(present(d))then call set_coe( dt, ext=d ) else call set_coe( dt, def=0.0 ) end if if(present(e))then call set_coe( et, ext=e ) else call set_coe( et, def=0.0 ) end if !-- 以下で先に格子間隔等の 1 回計算でよいものを求めておく. !-- これらは 1 方向のみで変化すればよい. !-- 格子点間隔の計算 do i=2,nx-1 dx(i)=(x(i+1)-x(i-1))*0.5 dx2(i)=dx(i)**2 end do do j=2,ny-1 dy(j)=(y(j+1)-y(j-1))*0.5 dy2(j)=dy(j)**2 end do dx(1)=(x(2)-x(1)) dx(nx)=(x(nx)-x(nx-1)) dy(1)=(y(2)-y(1)) dy(ny)=(y(ny)-y(ny-1)) do j=1,ny do i=1,nx dxdy(i,j)=dx(i)*dy(j) end do end do !-- ポアソン係数の計算 !-- 実際にポアソン係数が使用されるのは, 境界を除いた 1 回り内側なので, !-- 計算量削減のため, ループをそのようにしておく. !-- 最高次数係数 ac の計算 !$omp parallel default(shared) !$omp do private(i,j) do j=2,ny-1 do i=2,nx-1 ac(i,j)=2.0*(at(i,j)/dx2(i)+ct(i,j)/dy2(j)) adp(i,j)=(at(i,j)/(dx2(i))+0.5*dt(i,j)/dx(i))/ac(i,j) adm(i,j)=(at(i,j)/(dx2(i))-0.5*dt(i,j)/dx(i))/ac(i,j) cep(i,j)=(ct(i,j)/(dy2(j))+0.5*et(i,j)/dy(j))/ac(i,j) cem(i,j)=(ct(i,j)/(dy2(j))-0.5*et(i,j)/dy(j))/ac(i,j) bt(i,j)=0.25*bt(i,j)/(dxdy(i,j))/ac(i,j) end do end do !$omp end do !$omp end parallel err_max=eps ! while に入るための便宜的措置 !-- 実際のソルバ --- do while(err_max>=eps) err_max=0.0 do j=2,ny-1 do i=2,nx-1 tmp=-rho(i,j)/ac(i,j) divi(i,j)=1.0 !-- 以降, 反復計算に必要な点の周囲 8 点についてそれぞれ !-- inner_flag のチェックと同時に適切な値をそれぞれ逐次計算する, !-- 下で, 各方向の格子計算を行っているが, ib=8,-8 は 4 隅格子のときにしか !-- case select しない. なぜなら, 真横や上下に隅格子がくることはありえない. if(inner_flag(i,j).eqv..false.)then ! .false. なら領域計算開始 !-- 右横格子 if(inner_flag(i+1,j).eqv..false.)then tmp=tmp+adp(i,j)*psi(i+1,j) else select case (ib(i+1,j)) case(1) tmp=tmp+adp(i,j)*bnd(i+1,j) case(-2) ! 右横格子は 2 はありえない. tmp=tmp+adp(i,j)*bnd(i+1,j)*dx(i+1) divi(i,j)=divi(i,j)-adp(i,j) case(3) tmp=tmp+adp(i,j)*psi(2,j) case(4) tmp=tmp+adp(i,j)*(psi(i+1,j+1)-bnd(i+1,j)*dy(j)) case(-4) tmp=tmp+adp(i,j)*(psi(i+1,j-1)+bnd(i+1,j)*dy(j)) end select end if !-- 左横格子 if(inner_flag(i-1,j).eqv..false.)then tmp=tmp+adm(i,j)*psi(i-1,j) else select case (ib(i-1,j)) case(1) tmp=tmp+adm(i,j)*bnd(i-1,j) case(2) ! 左横格子は -2 はありえない. tmp=tmp-adm(i,j)*bnd(i-1,j)*dx(i-1) divi(i,j)=divi(i,j)-adm(i,j) case(3) tmp=tmp+adm(i,j)*psi(nx-1,j) case(4) tmp=tmp+adm(i,j)*(psi(i-1,j+1)-bnd(i-1,j)*dy(j)) case(-4) tmp=tmp+adm(i,j)*(psi(i-1,j-1)+bnd(i-1,j)*dy(j)) end select end if !-- 真上格子 if(inner_flag(i,j+1).eqv..false.)then tmp=tmp+cep(i,j)*psi(i,j+1) else select case (ib(i,j+1)) case(1) tmp=tmp+cep(i,j)*bnd(i,j+1) case(2) tmp=tmp+cep(i,j)*(psi(i+1,j+1)-bnd(i,j+1)*dx(i)) case(-2) tmp=tmp+cep(i,j)*(psi(i-1,j+1)+bnd(i,j+1)*dx(i)) case(3) tmp=tmp+cep(i,j)*psi(i,2) case(-4) ! 真上格子は 4 はありえない. tmp=tmp+cep(i,j)*bnd(i,j+1)*dy(j+1) divi(i,j)=divi(i,j)-cep(i,j) end select end if !-- 真下格子 if(inner_flag(i,j-1).eqv..false.)then tmp=tmp+cem(i,j)*psi(i,j-1) else select case (ib(i,j-1)) case(1) tmp=tmp+cem(i,j)*bnd(i,j-1) case(2) tmp=tmp+cem(i,j)*(psi(i+1,j-1)-bnd(i,j-1)*dx(i)) case(-2) tmp=tmp+cem(i,j)*(psi(i-1,j-1)+bnd(i,j-1)*dx(i)) case(3) tmp=tmp+cem(i,j)*psi(i,ny-1) case(4) ! 真下格子は -4 はありえない. tmp=tmp-cem(i,j)*bnd(i,j-1)*dy(j-1) divi(i,j)=divi(i,j)-cem(i,j) end select end if !-- 4 隅格子 if(signb==1)then ! そもそも bt = 0 なら計算しない. tmp_b=0.0 select case (ib(i-1,j-1)) ! 左下格子 case(1) tmp_b=tmp_b+bnd(i-1,j-1) case(2) tmp_b=tmp_b+psi(i,j-1)-bnd(i-1,j-1)*dx(i-1) case(-2) tmp_b=tmp_b+psi(i-2,j-1)+bnd(i-1,j-1)*dx(i-1) case(4) tmp_b=tmp_b+psi(i-1,j)-bnd(i-1,j-1)*dy(j-1) case(-4) tmp_b=tmp_b+psi(i-1,j-2)+bnd(i-1,j-1)*dy(j-1) case(3) ! 左下で周期境界なら, i==2 か j==2 しかありえない. if(i==2)then tmp_b=tmp_b+psi(nx-1,j-1) else if(j==2)then tmp_b=tmp_b+psi(i-1,ny-1) end if case(-3) ! 左下で -3 なら, 領域の左下隅しかない. tmp_b=tmp_b+psi(nx-1,ny-1) case(8) ! 左下では -8 は存在しない divi(i,j)=divi(i,j)-bt(i,j) tmp_b=tmp_b -0.5*(bnd(i,j-1)*dy(j-1)+bnd(i-1,j)*dx(i-1)) end select select case (ib(i+1,j+1)) ! 右上格子 case(1) tmp_b=tmp_b+bnd(i+1,j+1) case(2) tmp_b=tmp_b+psi(i,j+2)-bnd(i+1,j+1)*dx(i+1) case(-2) tmp_b=tmp_b+psi(i,j+1)+bnd(i+1,j+1)*dx(i+1) case(4) tmp_b=tmp_b+psi(i+1,j+2)-bnd(i+1,j+1)*dy(j+1) case(-4) tmp_b=tmp_b+psi(i+1,j)+bnd(i+1,j+1)*dy(j+1) case(3) ! 右上の場合, 周期境界は i==nx-1 か j==ny-1 しかない. if(i==nx-1)then tmp_b=tmp_b+psi(2,j-1) else if(j==ny-1)then tmp_b=tmp_b+psi(i-1,2) end if case(-3) ! 右上で -3 なら, 領域の右上隅しかない. tmp_b=tmp_b+psi(2,2) case(8) ! 右上では -8 は存在しない divi(i,j)=divi(i,j)-bt(i,j) tmp_b=tmp_b +0.5*(bnd(i,j+1)*dy(j+1)+bnd(i+1,j)*dx(i+1)) end select select case (ib(i-1,j+1)) ! 左上格子 case(1) tmp_b=tmp_b+bnd(i-1,j+1) case(2) tmp_b=tmp_b+psi(i,j+1)-bnd(i-1,j+1)*dx(i-1) case(-2) tmp_b=tmp_b+psi(i-2,j+1)+bnd(i-1,j+1)*dx(i-1) case(4) tmp_b=tmp_b+psi(i-1,j+2)-bnd(i-1,j+1)*dy(j+1) case(-4) tmp_b=tmp_b+psi(i-1,j)+bnd(i-1,j+1)*dy(j+1) case(3) ! 左上で周期境界なら, i==2 か j==ny-1 しかありえない. if(i==2)then tmp_b=tmp_b+psi(nx-1,j-1) else if(j==ny-1)then tmp_b=tmp_b+psi(i-1,2) end if case(-3) ! 左上で -3 なら, 領域の左上隅しかない. tmp_b=tmp_b+psi(nx-1,2) case(-8) ! 左上では 8 は存在しない divi(i,j)=divi(i,j)-bt(i,j) tmp_b=tmp_b +0.5*(bnd(i,j+1)*dy(j+1)-bnd(i-1,j)*dx(i-1)) end select select case (ib(i+1,j-1)) ! 右下格子 case(1) tmp_b=tmp_b+bnd(i+1,j-1) case(2) tmp_b=tmp_b+psi(i+2,j-1)-bnd(i+1,j-1)*dx(i+1) case(-2) tmp_b=tmp_b+psi(i,j-1)+bnd(i+1,j-1)*dx(i+1) case(4) tmp_b=tmp_b+psi(i-1,j)-bnd(i+1,j-1)*dy(j-1) case(-4) tmp_b=tmp_b+psi(i-1,j-2)+bnd(i+1,j-1)*dy(j-1) case(3) ! 右下で周期境界なら, i==nx-1 か j==2 しかありえない. if(i==nx-1)then tmp_b=tmp_b+psi(2,j-1) else if(j==2)then tmp_b=tmp_b+psi(i+1,ny-1) end if case(-3) ! 右下で -3 なら, 領域の右下隅しかない. tmp_b=tmp_b+psi(2,ny-1) case(-8) ! 右下では 8 は存在しない divi(i,j)=divi(i,j)-bt(i,j) tmp_b=tmp_b +0.5*(-bnd(i,j-1)*dy(j-1)+bnd(i+1,j)*dx(i+1)) end select end if tmp=tmp+bt(i,j)*tmp_b tmp=tmp/divi(i,j) end if err=abs(tmp-psi(i,j)) !-- 最大誤差の更新 if(err_max<=err)then err_max=err end if psi(i,j)=tmp end do end do end do !-- 境界の設定 call calculate_bound( ib, dx, dy, bnd, psi ) !-- 未定義領域には undef を代入する. do j=1,ny do i=1,nx if(ib(i,j)==10)then psi(i,j)=defun end if end do end do end subroutine Poisson_GauSei
Subroutine : | |||
x(:) : | real, intent(in)
| ||
y(:) : | real, intent(in)
| ||
rho(size(x),size(y)) : | real, intent(in)
| ||
eps : | real, intent(in)
| ||
boundary : | character(4), intent(in)
| ||
psi(size(x),size(y)) : | real, intent(inout)
| ||
bound_opt(size(x),size(y)) : | real, intent(in), optional
| ||
a(size(x),size(y)) : | real, intent(in), optional
| ||
b(size(x),size(y)) : | real, intent(in), optional
| ||
c(size(x),size(y)) : | real, intent(in), optional
| ||
d(size(x),size(y)) : | real, intent(in), optional
| ||
e(size(x),size(y)) : | real, intent(in), optional
| ||
undef : | real, intent(in), optional
| ||
inner_bound(size(x),size(y)) : | integer, intent(in), optional
| ||
init_flag : | logical, intent(in), optional
|
openmp によるスレッド並列が可能. ガウス・ザイデル法ではアルゴリズムの点から並列化が困難と思われたので, 並列計算によるポアソン方程式の求積が必要となるなら, ヤコビ法のものを使用されたい.
subroutine Poisson_Jacobi(x, y, rho, eps, boundary, psi, bound_opt, a, b, c, d, e, undef, inner_bound, init_flag ) ! openmp によるスレッド並列が可能. ! ガウス・ザイデル法ではアルゴリズムの点から並列化が困難と思われたので, ! 並列計算によるポアソン方程式の求積が必要となるなら, ! ヤコビ法のものを使用されたい. implicit none real, intent(in) :: x(:) ! 領域の横座標 real, intent(in) :: y(:) ! 領域の縦座標 real, intent(in) :: rho(size(x),size(y)) ! ポアソン方程式の強制項 ! rho =0 でラプラス方程式も求積可能 real, intent(in) :: eps ! 収束条件 character(4), intent(in) :: boundary ! 境界条件 ! 4 文字で各辺の境界条件を与える. ! 1 文字目 : x 下端, 2 文字目 : y 左端, 3 文字目 : x 上端, ! 4 文字目 : y 右端 ! boundary は 1 : 固定端境界, 2 : 自由端境界, 3 : 周期境界 real, intent(in), optional :: bound_opt(size(x),size(y)) ! 境界での強制 ! ノイマン境界の場合 : フラックス値 real, intent(in), optional :: a(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: b(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: c(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: d(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: e(size(x),size(y)) ! 各微分項の係数 real, intent(in), optional :: undef ! 未定義値 integer, intent(in), optional :: inner_bound(size(x),size(y)) ! 内部領域の境界. 値に応じてその格子点で境界値計算 ! 1 = 固定端境界, 10 = 境界の内側. ! 2 = y 方向自由端境界 (フラックスは上向き) ! -2 = y 方向自由端境界 (フラックスは下向き) ! 4 = x 方向自由端境界 (フラックスは右向き) ! -4 = x 方向自由端境界 (フラックスは左向き) ! 3 = 周期境界, -3 = 隅領域で両方とも周期境界 ! 8 = |_, ~| で両方とも自由境界条件 ! -8 = |~, _| で両方とも自由境界条件 ! この引数が与えられなければ全領域を計算する. ! 境界の内側格子点 (10) は反復計算を行わず, ! undef で設定された値もしくはゼロが入る. ! このときの境界値は bound_opt の値が用いられる. real, intent(inout) :: psi(size(x),size(y)) ! ポアソン方程式の解 logical, intent(in), optional :: init_flag ! psi の値をゼロで初期化するか. ! .true. = 初期化する. .false. = 初期化しない. ! デフォルトでは初期化する. integer :: i, j integer :: nx ! x 方向の配列要素 integer :: ny ! y 方向の配列要素 integer :: signb ! クロスターム b が存在するか integer, dimension(size(x),size(y)) :: ib real :: defun real :: err, err_max real :: bnd(size(x),size(y)) real :: dx(size(x)), dy(size(y)), dx2(size(x)), dy2(size(y)) real, dimension(size(x),size(y)) :: dxdy real, dimension(size(x),size(y)) :: at, bt, ct, dt, et real, dimension(size(x),size(y)) :: adp, adm, cep, cem, ac real, dimension(size(x),size(y)) :: tmp, tmp_b, divi character(4) :: bound logical, dimension(size(x),size(y)) :: inner_flag bound(1:4)=boundary(1:4) nx=size(x) ny=size(y) !-- 応答関数の初期化 if(present(init_flag))then if(init_flag.eqv..true.)then psi = 0.0 end if else psi = 0.0 end if !-- 内部境界の判別フラグの設定 if(present(inner_bound))then call set_bound( bound, ib, inner_flag, inner_bound ) else call set_bound( bound, ib, inner_flag ) end if !-- 領域・内部境界における境界値の設定 if(present(bound_opt))then call setval_bound( ib, bnd, psi, bound_opt ) else call setval_bound( ib, bnd, psi ) end if !-- 未定義値の設定 if(present(undef))then defun=undef else defun=0.0 end if !-- 係数の代入 !-- a, c については, 値が入れられていなければ, 全配列に 1 を代入する. if(present(a))then call set_coe( at, ext=a ) else call set_coe( at, def=1.0 ) end if if(present(c))then call set_coe( ct, ext=c ) else call set_coe( ct, def=1.0 ) end if if(present(b))then call set_coe( bt, ext=b ) signb=1 else call set_coe( bt, def=0.0 ) signb=0 end if if(present(d))then call set_coe( dt, ext=d ) else call set_coe( dt, def=0.0 ) end if if(present(e))then call set_coe( et, ext=e ) else call set_coe( et, def=0.0 ) end if !-- 以下で先に格子間隔等の 1 回計算でよいものを求めておく. !-- これらは 1 方向のみで変化すればよい. !-- 格子点間隔の計算 do i=2,nx-1 dx(i)=(x(i+1)-x(i-1))*0.5 dx2(i)=dx(i)**2 end do do j=2,ny-1 dy(j)=(y(j+1)-y(j-1))*0.5 dy2(j)=dy(j)**2 end do dx(1)=(x(2)-x(1)) dx(nx)=(x(nx)-x(nx-1)) dy(1)=(y(2)-y(1)) dy(ny)=(y(ny)-y(ny-1)) do j=1,ny do i=1,nx dxdy(i,j)=dx(i)*dy(j) end do end do !-- ポアソン係数の計算 !-- 実際にポアソン係数が使用されるのは, 境界を除いた 1 回り内側なので, !-- 計算量削減のため, ループをそのようにしておく. !-- 最高次数係数 ac の計算 !$omp parallel default(shared) !$omp do private(i,j) do j=2,ny-1 do i=2,nx-1 ac(i,j)=2.0*(at(i,j)/dx2(i)+ct(i,j)/dy2(j)) adp(i,j)=(at(i,j)/(dx2(i))+0.5*dt(i,j)/dx(i))/ac(i,j) adm(i,j)=(at(i,j)/(dx2(i))-0.5*dt(i,j)/dx(i))/ac(i,j) cep(i,j)=(ct(i,j)/(dy2(j))+0.5*et(i,j)/dy(j))/ac(i,j) cem(i,j)=(ct(i,j)/(dy2(j))-0.5*et(i,j)/dy(j))/ac(i,j) bt(i,j)=0.25*bt(i,j)/(dxdy(i,j))/ac(i,j) end do end do !$omp end do !$omp end parallel err_max=eps ! while に入るための便宜的措置 !-- 実際のソルバ --- do while(err_max>=eps) err_max=0.0 !$omp parallel default(shared) !$omp do schedule(dynamic) private(i,j) do j=2,ny-1 do i=2,nx-1 tmp(i,j)=-rho(i,j)/ac(i,j) divi(i,j)=1.0 !-- 以降, 反復計算に必要な点の周囲 8 点についてそれぞれ !-- inner_flag のチェックと同時に適切な値をそれぞれ逐次計算する, !-- 下で, 各方向の格子計算を行っているが, ib=8,-8 は 4 隅格子のときにしか !-- case select しない. なぜなら, 真横や上下に隅格子がくることはありえない. if(inner_flag(i,j).eqv..false.)then ! .false. なら領域計算開始 !-- 右横格子 if(inner_flag(i+1,j).eqv..false.)then tmp(i,j)=tmp(i,j)+adp(i,j)*psi(i+1,j) else select case (ib(i+1,j)) case(1) tmp(i,j)=tmp(i,j)+adp(i,j)*bnd(i+1,j) case(-2) ! 右横格子は 2 はありえない. tmp(i,j)=tmp(i,j)+adp(i,j)*bnd(i+1,j)*dx(i+1) divi(i,j)=divi(i,j)-adp(i,j) case(3) tmp(i,j)=tmp(i,j)+adp(i,j)*psi(2,j) case(4) tmp(i,j)=tmp(i,j)+adp(i,j)*(psi(i+1,j+1)-bnd(i+1,j)*dy(j)) case(-4) tmp(i,j)=tmp(i,j)+adp(i,j)*(psi(i+1,j-1)+bnd(i+1,j)*dy(j)) end select end if !-- 左横格子 if(inner_flag(i-1,j).eqv..false.)then tmp(i,j)=tmp(i,j)+adm(i,j)*psi(i-1,j) else select case (ib(i-1,j)) case(1) tmp(i,j)=tmp(i,j)+adm(i,j)*bnd(i-1,j) case(2) ! 左横格子は -2 はありえない. tmp(i,j)=tmp(i,j)-adm(i,j)*bnd(i-1,j)*dx(i-1) divi(i,j)=divi(i,j)-adm(i,j) case(3) tmp(i,j)=tmp(i,j)+adm(i,j)*psi(nx-1,j) case(4) tmp(i,j)=tmp(i,j)+adm(i,j)*(psi(i-1,j+1)-bnd(i-1,j)*dy(j)) case(-4) tmp(i,j)=tmp(i,j)+adm(i,j)*(psi(i-1,j-1)+bnd(i-1,j)*dy(j)) end select end if !-- 真上格子 if(inner_flag(i,j+1).eqv..false.)then tmp(i,j)=tmp(i,j)+cep(i,j)*psi(i,j+1) else select case (ib(i,j+1)) case(1) tmp(i,j)=tmp(i,j)+cep(i,j)*bnd(i,j+1) case(2) tmp(i,j)=tmp(i,j)+cep(i,j)*(psi(i+1,j+1)-bnd(i,j+1)*dx(i)) case(-2) tmp(i,j)=tmp(i,j)+cep(i,j)*(psi(i-1,j+1)+bnd(i,j+1)*dx(i)) case(3) tmp(i,j)=tmp(i,j)+cep(i,j)*psi(i,2) case(-4) ! 真上格子は 4 はありえない. tmp(i,j)=tmp(i,j)+cep(i,j)*bnd(i,j+1)*dy(j+1) divi(i,j)=divi(i,j)-cep(i,j) end select end if !-- 真下格子 if(inner_flag(i,j-1).eqv..false.)then tmp(i,j)=tmp(i,j)+cem(i,j)*psi(i,j-1) else select case (ib(i,j-1)) case(1) tmp(i,j)=tmp(i,j)+cem(i,j)*bnd(i,j-1) case(2) tmp(i,j)=tmp(i,j)+cem(i,j)*(psi(i+1,j-1)-bnd(i,j-1)*dx(i)) case(-2) tmp(i,j)=tmp(i,j)+cem(i,j)*(psi(i-1,j-1)+bnd(i,j-1)*dx(i)) case(3) tmp(i,j)=tmp(i,j)+cem(i,j)*psi(i,ny-1) case(4) ! 真下格子は -4 はありえない. tmp(i,j)=tmp(i,j)-cem(i,j)*bnd(i,j-1)*dy(j-1) divi(i,j)=divi(i,j)-cem(i,j) end select end if !-- 4 隅格子 if(signb==1)then ! そもそも bt = 0 なら計算しない. tmp_b(i,j)=0.0 select case (ib(i-1,j-1)) ! 左下格子 case(1) tmp_b(i,j)=tmp_b(i,j)+bnd(i-1,j-1) case(2) tmp_b(i,j)=tmp_b(i,j)+psi(i,j-1)-bnd(i-1,j-1)*dx(i-1) case(-2) tmp_b(i,j)=tmp_b(i,j)+psi(i-2,j-1)+bnd(i-1,j-1)*dx(i-1) case(4) tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j)-bnd(i-1,j-1)*dy(j-1) case(-4) tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j-2)+bnd(i-1,j-1)*dy(j-1) case(3) ! 左下で周期境界なら, i==2 か j==2 しかありえない. if(i==2)then tmp_b(i,j)=tmp_b(i,j)+psi(nx-1,j-1) else if(j==2)then tmp_b(i,j)=tmp_b(i,j)+psi(i-1,ny-1) end if case(-3) ! 左下で -3 なら, 領域の左下隅しかない. tmp_b(i,j)=tmp_b(i,j)+psi(nx-1,ny-1) case(8) ! 左下では -8 は存在しない divi(i,j)=divi(i,j)-bt(i,j) tmp_b(i,j)=tmp_b(i,j) -0.5*(bnd(i,j-1)*dy(j-1)+bnd(i-1,j)*dx(i-1)) end select select case (ib(i+1,j+1)) ! 右上格子 case(1) tmp_b(i,j)=tmp_b(i,j)+bnd(i+1,j+1) case(2) tmp_b(i,j)=tmp_b(i,j)+psi(i,j+2)-bnd(i+1,j+1)*dx(i+1) case(-2) tmp_b(i,j)=tmp_b(i,j)+psi(i,j+1)+bnd(i+1,j+1)*dx(i+1) case(4) tmp_b(i,j)=tmp_b(i,j)+psi(i+1,j+2)-bnd(i+1,j+1)*dy(j+1) case(-4) tmp_b(i,j)=tmp_b(i,j)+psi(i+1,j)+bnd(i+1,j+1)*dy(j+1) case(3) ! 右上の場合, 周期境界は i==nx-1 か j==ny-1 しかない. if(i==nx-1)then tmp_b(i,j)=tmp_b(i,j)+psi(2,j-1) else if(j==ny-1)then tmp_b(i,j)=tmp_b(i,j)+psi(i-1,2) end if case(-3) ! 右上で -3 なら, 領域の右上隅しかない. tmp_b(i,j)=tmp_b(i,j)+psi(2,2) case(8) ! 右上では -8 は存在しない divi(i,j)=divi(i,j)-bt(i,j) tmp_b(i,j)=tmp_b(i,j) +0.5*(bnd(i,j+1)*dy(j+1)+bnd(i+1,j)*dx(i+1)) end select select case (ib(i-1,j+1)) ! 左上格子 case(1) tmp_b(i,j)=tmp_b(i,j)+bnd(i-1,j+1) case(2) tmp_b(i,j)=tmp_b(i,j)+psi(i,j+1)-bnd(i-1,j+1)*dx(i-1) case(-2) tmp_b(i,j)=tmp_b(i,j)+psi(i-2,j+1)+bnd(i-1,j+1)*dx(i-1) case(4) tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j+2)-bnd(i-1,j+1)*dy(j+1) case(-4) tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j)+bnd(i-1,j+1)*dy(j+1) case(3) ! 左上で周期境界なら, i==2 か j==ny-1 しかありえない. if(i==2)then tmp_b(i,j)=tmp_b(i,j)+psi(nx-1,j-1) else if(j==ny-1)then tmp_b(i,j)=tmp_b(i,j)+psi(i-1,2) end if case(-3) ! 左上で -3 なら, 領域の左上隅しかない. tmp_b(i,j)=tmp_b(i,j)+psi(nx-1,2) case(-8) ! 左上では 8 は存在しない divi(i,j)=divi(i,j)-bt(i,j) tmp_b(i,j)=tmp_b(i,j) +0.5*(bnd(i,j+1)*dy(j+1)-bnd(i-1,j)*dx(i-1)) end select select case (ib(i+1,j-1)) ! 右下格子 case(1) tmp_b(i,j)=tmp_b(i,j)+bnd(i+1,j-1) case(2) tmp_b(i,j)=tmp_b(i,j)+psi(i+2,j-1)-bnd(i+1,j-1)*dx(i+1) case(-2) tmp_b(i,j)=tmp_b(i,j)+psi(i,j-1)+bnd(i+1,j-1)*dx(i+1) case(4) tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j)-bnd(i+1,j-1)*dy(j-1) case(-4) tmp_b(i,j)=tmp_b(i,j)+psi(i-1,j-2)+bnd(i+1,j-1)*dy(j-1) case(3) ! 右下で周期境界なら, i==nx-1 か j==2 しかありえない. if(i==nx-1)then tmp_b(i,j)=tmp_b(i,j)+psi(2,j-1) else if(j==2)then tmp_b(i,j)=tmp_b(i,j)+psi(i+1,ny-1) end if case(-3) ! 右下で -3 なら, 領域の右下隅しかない. tmp_b(i,j)=tmp_b(i,j)+psi(2,ny-1) case(-8) ! 右下では 8 は存在しない divi(i,j)=divi(i,j)-bt(i,j) tmp_b(i,j)=tmp_b(i,j) +0.5*(-bnd(i,j-1)*dy(j-1)+bnd(i+1,j)*dx(i+1)) end select end if tmp(i,j)=tmp(i,j)+bt(i,j)*tmp_b(i,j) tmp(i,j)=tmp(i,j)/divi(i,j) end if end do end do !$omp end do !$omp end parallel !!-- ここまでは if 文のインデントを調節しない. ! end if ! ! tmp=tmp/ac(i,j) !-- 誤差の計算 --- do j=2,ny-1 do i=2,nx-1 err=abs(tmp(i,j)-psi(i,j)) !-- 最大誤差の更新 if(err_max<=err)then err_max=err end if end do end do !-- 一斉更新 do j=2,ny-1 do i=2,nx-1 psi(i,j)=tmp(i,j) end do end do end do !-- 境界の設定 call calculate_bound( ib, dx, dy, bnd, psi ) !-- 未定義領域には undef を代入する. do j=1,ny do i=1,nx if(ib(i,j)==10)then psi(i,j)=defun end if end do end do end subroutine Poisson_Jacobi
Subroutine : | |||
ib(:,:) : | integer, intent(in)
| ||
dx(size(ib,1)) : | real, intent(in)
| ||
dy(size(ib,2)) : | real, intent(in)
| ||
bnd(size(ib,1),size(ib,2)) : | real, intent(in)
| ||
psi(size(ib,1),size(ib,2)) : | real, intent(inout)
|
ib 判別を元に, ノイマン型, 周期境界型について境界値を計算する.
subroutine calculate_bound( ib, dx, dy, bnd, psi ) ! ib 判別を元に, ノイマン型, 周期境界型について境界値を計算する. integer, intent(in) :: ib(:,:) ! 境界整数判別 real, intent(in) :: dx(size(ib,1)) ! x 方向の格子解像度 real, intent(in) :: dy(size(ib,2)) ! y 方向の格子解像度 real, intent(in) :: bnd(size(ib,1),size(ib,2)) ! 境界値(ノイマン型のみ使用) real, intent(inout) :: psi(size(ib,1),size(ib,2)) ! 応答関数 integer :: i, j, nx, ny nx=size(ib,1) ny=size(ib,2) do j=1,ny do i=1,nx select case (ib(i,j)) case (1) psi(i,j)=bnd(i,j) case (2) ! x 方向にフラックス一定, 上側が参照値 psi(i,j)=psi(i+1,j)-bnd(i,j)*dx(i) case (-2) ! x 方向にフラックス一定, 下側が参照値 psi(i,j)=psi(i-1,j)+bnd(i,j)*dx(i) case (4) ! y 方向にフラックス一定, 右側が参照値 psi(i,j)=psi(i,j+1)-bnd(i,j)*dy(j) case (-4) ! y 方向にフラックス一定, 左側が参照値 psi(i,j)=psi(i,j-1)+bnd(i,j)*dy(j) case (3) ! 周期境界 if(j==1)then psi(i,j)=psi(i,ny-1) else if(j==ny)then psi(i,j)=psi(i,2) else if(i==1)then psi(i,j)=psi(nx-1,j) else if(i==nx)then psi(i,j)=psi(2,j) end if case (-3) ! 4 隅で同時周期境界 if(i==1.and.j==1)then psi(i,j)=psi(nx-1,ny-1) else if(i==1.and.j==ny)then psi(i,j)=psi(nx-1,2) else if(i==nx.and.j==1)then psi(i,j)=psi(2,ny-1) else if(i==nx.and.j==ny)then psi(i,j)=psi(2,2) end if case (8) ! 両方フラックス一定で左下角か右上角 if(i==1.and.j==1)then ! -- 評価 1 psi(i,j)=psi(i+1,j+1)-0.5*(bnd(i+1,j)*dy(j)+bnd(i,j+1)*dx(i)) else if(i==nx.and.j==ny)then ! -- 評価 2 psi(i,j)=psi(i-1,j-1)+0.5*(bnd(i-1,j)*dy(j)+bnd(i,j-1)*dx(i)) else if(ib(i-1,j)==10.or.ib(i,j-1)==10)then ! -- 評価 1 と同じ psi(i,j)=psi(i+1,j+1)-0.5*(bnd(i+1,j)*dy(j)+bnd(i,j+1)*dx(i)) else if(ib(i+1,j)==10.or.ib(i,j+1)==10)then ! -- 評価 2 と同じ psi(i,j)=psi(i-1,j-1)+0.5*(bnd(i-1,j)*dy(j)+bnd(i,j-1)*dx(i)) end if case (-8) ! 両方フラックス一定で右下角か左上角 if(i==1.and.j==ny)then ! -- 評価 1 psi(i,j)=psi(i+1,j+1)+0.5*(bnd(i+1,j)*dy(j)-bnd(i,j-1)*dx(i)) else if(i==nx.and.j==1)then ! -- 評価 2 psi(i,j)=psi(i-1,j-1)+0.5*(-bnd(i-1,j)*dy(j)+bnd(i,j+1)*dx(i)) else if(ib(i-1,j)==10.or.ib(i,j+1)==10)then ! -- 評価 1 と同じ psi(i,j)=psi(i+1,j+1)+0.5*(bnd(i+1,j)*dy(j)-bnd(i,j-1)*dx(i)) else if(ib(i+1,j)==10.or.ib(i,j-1)==10)then ! -- 評価 2 と同じ psi(i,j)=psi(i-1,j-1)+0.5*(-bnd(i-1,j)*dy(j)+bnd(i,j+1)*dx(i)) end if end select end do end do end subroutine calculate_bound
Subroutine : | |||
bound : | character(4), intent(in)
| ||
ib(:,:) : | integer, intent(inout)
| ||
inner_flag(size(ib,1),size(ib,2)) : | logical, intent(inout)
| ||
inner_bound(size(ib,1),size(ib,2)) : | integer, intent(in), optional
|
各反復法ルーチンにおいて設定される境界条件のフラグをチェック, 設定する.
subroutine set_bound( bound, ib, inner_flag, inner_bound ) ! 各反復法ルーチンにおいて設定される境界条件のフラグをチェック, 設定する. implicit none character(4), intent(in) :: bound ! 領域境界のフラグ integer, intent(inout) :: ib(:,:) ! 全境界の判別整数 logical, intent(inout) :: inner_flag(size(ib,1),size(ib,2)) ! 全領域境界フラグ integer, intent(in), optional :: inner_bound(size(ib,1),size(ib,2)) ! 内部領域境界の判別整数 integer :: i, j, nx, ny nx=size(ib,1) ny=size(ib,2) !-- 出力変数の初期化 ib=0 inner_flag=.false. !-- 周期境界の設定確認. !-- 周期境界なので, 両端とも 3 が設定されていないといけない. if(bound(1:1)=='3')then if(bound(3:3)/='3')then write(*,*) "### ERROR ###" write(*,*) "if bound = 3, bound(1:1)==bound(3:3). STOP." stop end if end if if(bound(3:3)=='3')then if(bound(1:1)/='3')then write(*,*) "### ERROR ###" write(*,*) "if bound = 3, bound(1:1)==bound(3:3). STOP." stop end if end if if(bound(2:2)=='3')then if(bound(4:4)/='3')then write(*,*) "### ERROR ###" write(*,*) "if bound = 3, bound(2:2)==bound(4:4). STOP." stop end if end if if(bound(4:4)=='3')then if(bound(2:2)/='3')then write(*,*) "### ERROR ###" write(*,*) "if bound = 3, bound(2:2)==bound(4:4). STOP." stop end if end if select case (bound(1:1)) case ('1') do i=2,nx-1 ib(i,1)=1 end do case ('2') do i=2,nx-1 ib(i,1)=4 ! y 方向のフラックスで参照値は上側 end do case ('3') do i=2,nx-1 ib(i,1)=3 ! y 方向 (x 軸) 周期境界 end do end select select case (bound(2:2)) case ('1') do j=2,ny-1 ib(1,j)=1 end do case ('2') do j=2,ny-1 ib(1,j)=2 ! x 方向のフラックスで参照値は右側 end do case ('3') do j=2,ny-1 ib(1,j)=-3 ! x 方向 (y 軸) 周期境界 end do end select select case (bound(3:3)) case ('1') do i=2,nx-1 ib(i,ny)=1 end do case ('2') do i=2,nx-1 ib(i,ny)=-4 ! y 方向のフラックスで参照値は下側 end do case ('3') do i=2,nx-1 ib(i,ny)=3 ! y 方向 (x 軸) 周期境界 end do end select select case (bound(4:4)) case ('1') do j=2,ny-1 ib(nx,j)=1 end do case ('2') do j=2,ny-1 ib(nx,j)=-2 ! x 方向のフラックスで参照値は左側 end do case ('3') do j=2,ny-1 ib(nx,j)=-3 ! x 方向 (y 軸) 周期境界 end do end select !-- 領域隅境界は 2 辺が重なるので, 境界での強制方法は以下の順に優先する. !-- (1) どちらかが 1 のとき, 隅領域は 1, !-- (2) (1) 以外でどちらかが 3 のとき, 隅領域は 3. !-- (3) (1), (2) 以外で 2. つまり, 両軸とも 2 なら 2 となる. ib(1,1)=0 ib(1,ny)=0 ib(nx,1)=0 ib(nx,ny)=0 if(bound(1:1)=='1')then ib(1,1)=1 ib(nx,1)=1 end if if(bound(2:2)=='1')then ib(1,1)=1 ib(1,ny)=1 end if if(bound(3:3)=='1')then ib(1,ny)=1 ib(nx,ny)=1 end if if(bound(4:4)=='1')then ib(nx,1)=1 ib(nx,ny)=1 end if !-- 4 隅とも周期境界の場合 if(bound(1:2)=='33')then ib(1,1)=-3 ib(1,ny)=-3 ib(nx,1)=-3 ib(nx,ny)=-3 end if if(bound(1:1)=='3')then if(ib(1,1)==0)then ib(1,1)=3 end if if(ib(nx,1)==0)then ib(nx,1)=3 end if end if if(bound(2:2)=='3')then if(ib(1,1)==0)then ib(1,1)=3 end if if(ib(1,ny)==0)then ib(1,ny)=3 end if end if if(bound(3:3)=='3')then if(ib(1,ny)==0)then ib(1,ny)=3 end if if(ib(nx,ny)==0)then ib(nx,ny)=3 end if end if if(bound(4:4)=='3')then if(ib(nx,1)==0)then ib(nx,1)=3 end if if(ib(nx,ny)==0)then ib(nx,ny)=3 end if end if !-- 領域の隅であることを 8, -8 で示す. !-- 8 の場合, 左下か右上であることを, -8 の場合, 左上か右下であることを示す. if(ib(1,1)==0)then ib(1,1)=8 end if if(ib(nx,1)==0)then ib(nx,1)=-8 end if if(ib(1,ny)==0)then ib(1,ny)=-8 end if if(ib(nx,ny)==0)then ib(nx,ny)=8 end if !-- 内部境界の設定 if(present(inner_bound))then do j=2,ny-1 do i=2,nx-1 ib(i,j)=inner_bound(i,j) if(inner_bound(i,j)/=0)then inner_flag(i,j)=.true. end if end do end do end if do i=1,nx ! inner_bound が設定されているに関わらず, 領域端は計算しない. inner_flag(i,1)=.true. inner_flag(i,ny)=.true. end do do j=1,ny inner_flag(1,j)=.true. inner_flag(nx,j)=.true. end do end subroutine set_bound
Subroutine : | |||
coe(:,:) : | real, intent(inout)
| ||
ext(size(coe,1),size(coe,2)) : | real, intent(in), optional
| ||
def : | real, intent(in), optional
|
2 次元配列に ext で指定された値もしくは def で指定された一定値を代入する. ext, def どちらも optional であるが, 必ずどちらかは指定されていないといけない.
subroutine set_coe( coe, ext, def ) ! 2 次元配列に ext で指定された値もしくは def で指定された一定値を代入する. ! ext, def どちらも optional であるが, 必ずどちらかは指定されていないといけない. implicit none real, intent(inout) :: coe(:,:) ! 代入される配列 real, intent(in), optional :: ext(size(coe,1),size(coe,2)) ! 代入する配列 real, intent(in), optional :: def ! 代入する一定値 integer :: i, j, nx, ny nx=size(coe,1) ny=size(coe,2) if(present(ext))then do j=1,ny do i=1,nx coe(i,j)=ext(i,j) end do end do else if(present(def))then do j=1,ny do i=1,nx coe(i,j)=def end do end do else write(*,*) "### ERROR ###" write(*,*) "subroutine set_coe must be set optional argument 'ext' or 'def'" write(*,*) "STOP." stop end if end subroutine set_coe
Subroutine : | |||
ib(:,:) : | integer, intent(in)
| ||
bnd(size(ib,1),size(ib,2)) : | real, intent(inout)
| ||
psi(size(ib,1),size(ib,2)) : | real, intent(inout)
| ||
bound_opt(size(ib,1),size(ib,2)) : | real, intent(in), optional
|
境界値を境界条件判別整数をもとに設定する.
subroutine setval_bound( ib, bnd, psi, bound_opt ) ! 境界値を境界条件判別整数をもとに設定する. implicit none integer, intent(in) :: ib(:,:) ! 境界条件判別整数 real, intent(inout) :: bnd(size(ib,1),size(ib,2)) ! 境界での値 real, intent(inout) :: psi(size(ib,1),size(ib,2)) ! 応答 real, intent(in), optional :: bound_opt(size(ib,1),size(ib,2)) ! 境界での値 integer :: i, j, nx, ny nx=size(ib,1) ny=size(ib,2) !-- 境界値の設定 if(present(bound_opt))then do j=1,ny do i=1,nx bnd(i,j)=bound_opt(i,j) end do end do else do j=1,ny do i=1,nx bnd(i,j)=0.0 end do end do end if !-- 境界条件の代入 "ib(i,j)==1 の場合のみ" !-- 内部領域についてもここで代入してしまう. do j=1,ny do i=1,nx if(ib(i,j)==1)then psi(i,j)=bnd(i,j) end if end do end do end subroutine setval_bound