@zigen 's note

CIP Memo

最終更新:

mynote

- view
管理者のみ編集可

CIP とは

スプライン近似に似た(3次元で補間する)もの
非移流項が入ってくる場合はSemi-Lagrangeで解く

ラグランジュ補間とスプライン補間の違いについては下のHPを参照


CIP法とは
東京工業大学の矢部孝教授1980年代に考案した数値計算手法。
Particle In Cell Scheme(PIC法)の研究を行っていた
↓・・・が嫌気がさし、自分が開発したスキームに皮肉を込めてPICを逆から読んで。
Cubic-Interpolated Pseudo-Particle Schemeと名付けた。CIPとはその略である。
(Cubic-Interpolatedの由来は四で字のごとく「3次元で補間」の意味)
↓しかし、三次元補間なら他にも在るわけで・・・
Constrained Interpolation Profile Schemeに名前変更
(この辺は完全に言葉遊びらしいです(この前講演してくださったときに言ってました(笑)))

  • よって論文検索では「CIP」で検索してもいいですが、正式名称は「Cubic-Interpolated Pseudo-Particle Scheme」「Constrained Interpolation Profile Scheme」の二種類在ります!!
    • リスト1
    • リスト2

テストコード

program main
implicit none
double precision :: dt, dx, o_dx, o_dx2, v_v, f_d, dx_isign, xx, a, b
double precision :: dft, dftx
double precision, dimension(:), allocatable :: f, fn, g, gn, ft, ftx, initial
integer :: i, j, t, x_up, n_t = 3, n_x = 20

allocate( initial(0:n_x+1), f(0:n_x+1), fn(0:n_x+1), g(0:n_x+1), gn(0:n_x+1), ft(0:n_x+1), ftx(0:n_x+1) )

!###########< Set initial condition (start) >##########
!
  dx = 1.d0/dble(n_x)
  o_dx = 1.d0/dx
  o_dx2 = o_dx*o_dx
  dt = 1.d0/n_t
  dx_isign = 1.d0
  v_v = 0.1d0
  xx = -v_v*dt
  dx_isign = -dsign(1.d0, v_v)

  do i = 0, n_x+1
     initial(i) = 0.d0
     f(i)  = 0.d0
     fn(i) = 0.d0
     ft(i) = 0.d0
     ftx(i) = 0.d0
     gn(i) = 0.d0
        if( 5<=i .and. i<=10 ) f(i) = 1.d0
        if( 5<=i .and. i<=10 ) initial(i) = 1.d0
     enddo
     do i = 0, n_x
        g(i) = (f(i) + f(i+1))*o_dx
     enddo
        g(n_x+1) = 0.d0
!
!###########< Set initial condition (end) >##########
!=======================================================================
! Main Loop
!=======================================================================
  do t = 1, n_t
     do i = 1, n_x
        x_up = i + int(dx_isign)
        f_d = (f(x_up) - f(i))*o_dx*dx_isign
        a = (g(i) + g(x_up) - 2.d0*f_d)*o_dx2
        b = (3.d0*f_d - 2.d0*g(i) - g(x_up))*o_dx*dx_isign
        fn(i) = ((a*xx + b)*xx + g(i))*xx + f(i)
        gn(i) = (3.d0*a*xx + 2.d0*b)*xx + g(i)
     enddo
        f(1:n_x) = fn(1:n_x)
        g(1:n_x) = gn(1:n_x)
        f(0) = fn(n_x+1)
        g(0) = gn(n_x+1)
  enddo

  open(99, file="result_cip.txt")
     do i = 1, n_x
        write(99,*) i, initial(i), f(i)
        write(*,*) i, initial(i), f(i)
     enddo
  close(99)
end program main
目安箱バナー