fortran66のブログ

fortran について書きます。

コンデンサの電位を計算

パリティ Vol.06 No.03 (1991) 62. のフォートラン入門より、平行平板コンデンサの作る二次元電場の電位を、ラプラス方程式のヤコビによる解法で解きます。

forallを用い mask 配列を使ってコンデンサの電極は計算を飛ばすことにします。forall は並列代入操作で、ただの do loop の時と配列更新のタイミングが違い、値は全ループが終わるまで更新されません。

さぼって収束性は確かめずに、十分と思われる回数反復することにします。計算の核になる部分はほんとに短く書けます。

ラプラス方程式は、熱伝導やゴムを伸ばした時の形状も記述しているのでコンデンサの電極の代わりに自由に定値を与えてやれば、ラバーフェチも満足の計算ができます。

グラフィックは、OpenGLイタレーションごとに描きたかったのですが、めんどくさかったので収束後の電位だけを、いい加減にプロットすることにしましたw

実行結果

ソース・プログラム

    program jacobi
      use plotter
      implicit none
      integer, parameter :: nsize = 50, nsz = nsize + 1
      real :: v(-nsz:nsz, -nsz:nsz) = 0.0
      logical :: mask(-nsz:nsz, -nsz:nsz) = .true.  
      integer :: i, j, iter
      mask(-10, -5:5) = .false. ! static voltage
      mask( 10, -5:5) = .false. 
      v(-10, -5:5) =  10.0 
      v( 10, -5:5) = -10.0
      !
      do iter = 1, 200 
        forall(i = -nsize:nsize, j = -nsize:nsize, mask(i, j)) &
             v(i, j) = 0.25 * ( v(i - 1, j) + v(i + 1, j) + v(i, j - 1) + v(i, j + 1) )
      end do   
      ! 
      call graph()
      stop
    
    contains 
      
      subroutine graph()
        integer :: ix, iy
        real :: x, y
        call gr_on('Laplace equation', 600, 600)
        call gr_pencol(0)
        ! x-direction     
        do iy = -nsize, nsize, 2
          x = 0.5 + 0.005 * -nsize + 0.005 * iy
          y = 0.5 + 0.003 * iy
          call move( x, y )
          do ix = -nsize, nsize, 2
            x = 0.005 * ix + 0.005 * iy + 0.5 
            y = 0.5 + 0.003 * iy + v(ix, iy) * 0.02
            call draw( x, y )
          end do 
        end do
        ! y-direction
        do ix = -nsize, nsize, 2 
          x = 0.005 * ix + 0.005 * -nsize + 0.5 
          y = 0.5 + 0.003 * -nsize
          call move( x, y )
          do iy = -nsize, nsize, 2
            x = 0.005 * ix + 0.005 * iy + 0.5 
            y = 0.5 + 0.003 * iy + v(ix, iy) * 0.02
            call draw( x, y )
          end do 
        end do
        call gr_off()
        return
      end subroutine graph
      
    end program jacobi

グラフィックス・ルーチン

module uho_win
  use ifwina
  use ifwinty
  use ifmt, only : RTL_CRITICAL_SECTION
  implicit none
  type :: t_wnd
    integer (HANDLE) :: hWnd
    integer (HANDLE) :: hDC
    integer (LPINT)  :: hThread
    integer         :: id
    integer (HANDLE) :: hPen
    character (LEN = 80) :: title = 'Fortran Plot'
    integer         :: nsize_x = 320, nsize_y = 240  
    real            :: xmin = 0.0, xmax = 1.0, ymin = 0.0, ymax = 1.0   
  end type
  type (t_wnd) :: wnd
  type (RTL_CRITICAL_SECTION) :: lpCriticalSection
 contains
 !--------------------------------------------------------------------------------
  integer(4) function WinMain( hInstance, hPrevInstance, lpszCmdLine, nCmdShow, wnd)
    implicit none
    integer (HANDLE), intent(in) :: hInstance, hPrevInstance 
    integer (LPSTR) , intent(in) :: lpszCmdLine
    integer (SINT)  , intent(in) :: nCmdShow
    type (t_wnd)    , intent(in out) :: wnd
    ! Variables
    logical, save :: qfirst = .true.
    type (T_WNDCLASS) :: wc
    type (T_MSG)      :: mesg
    integer (HANDLE)  :: hWndMain
    integer (BOOL)    :: iretb
    character (LEN = 256) :: ClassName = 'Fortran'//char(0)
    integer :: noffset_x, noffset_y !side line = 6, title bar = 25
    !  Init   Main window
    noffset_x = 2 * GetSystemMetrics(SM_CXFIXEDFRAME)
    noffset_y = 2 * GetSystemMetrics(SM_CYFIXEDFRAME) + GetSystemMetrics(SM_CYCAPTION)
    !
    WinMain = -1 ! Error code 
    if (qfirst) then
      wc%lpszClassName =  loc(ClassName)
      wc%lpfnWndProc   =  loc(MainWndProc)               ! call BACK procedure name
      wc%style        = ior(CS_VREDRAW , CS_HREDRAW)
      wc%hInstance     = hInstance
      wc%hIcon        = NULL   
      wc%hCursor      = LoadCursor( NULL, IDC_ARROW )
      wc%hbrBackground = ( COLOR_WINDOW + 1 )
      if ( RegisterClass(wc) == 0 ) return
      qfirst = .FALSE.
    end if
    !Init instance
    WinMain = -2 ! Error code 
    hWndMain = CreateWindowEx(  0, ClassName,                        &
                                trim(wnd%title)//char(0),           &
                                int(ior(WS_OVERLAPPED, WS_SYSMENU)), &
                                CW_USEDEFAULT, CW_USEDEFAULT,        &
                                wnd%nsize_x + noffset_x,            &
                                wnd%nsize_y + noffset_y,            &
                                0,                                &
                                0,                                &
                                hInstance,                         &
                                 loc(wnd)                           ) 
    if (hWndMain == 0) return
    iretb = ShowWindow( hWndMain, nCmdShow )
    iretb = UpdateWindow( hWndMain )
    ! Message Loop
    do while ( GetMessage (mesg, NULL, 0, 0) ) 
      iretb = TranslateMessage( mesg ) 
      iretb = DispatchMessage(  mesg )
    end do
    WinMain = mesg%wParam
    return
  end function WinMain
  !----------------------------------------------------------
  integer (LRESULT) function MainWndProc( hWnd, mesg, wParam, lParam ) 
  !DEC$ ATTRIBUTES STDcall, DECORATE, ALIAS : 'MainWndProc' :: MainWndProc
    use IFMT
    implicit none
    integer (HANDLE) , intent(in) :: hWnd
    integer (UINT)   , intent(in) :: mesg
    integer (fwParam), intent(in) :: wParam
    integer (flParam), intent(in) :: lParam
    !
    integer (HANDLE) :: hDC, hBmp
    integer (BOOL)   :: iretb
    type (T_PAINTSTRUCT)   :: ps
    type (T_RECT)         :: rc
    !
    MainWndProc = 0
    select case ( mesg )
      case (WM_CREATE)
       hDC      = GetDC(hWnd)
       iretb    = GetClientRect(hWnd, rc)
       hBmp     = CreateCompatibleBitmap(hDC, rc%right - rc%left, rc%bottom - rc%top)
       wnd%hWnd = hWnd
       wnd%hDC  = CreateCompatibleDC(hDC)
       iretb    = SelectObject(wnd%hDC, hBmp)
       iretb    = PatBlt(wnd%hDC, 0, 0, rc%right - rc%left, rc%bottom - rc%top, WHITENESS)
       iretb    = ReleaseDC(hWnd, hDC)
       iretb    = DeleteObject(hBmp)
     case (WM_DESTROY)
       iretb = DeleteObject( wnd%hDC )
       call PostQuitMessage( 0 )
     case (WM_PAINT)
       call EnterCriticalSection(  loc(lpCriticalSection) )
       hDC    = BeginPaint(    wnd%hWnd, ps )
       iretb  = GetClientRect( wnd%hWnd, rc )
       iretb  = BitBlt(hDC, 0, 0, rc%right  - rc%left, rc%bottom - rc%top, wnd%hDC, 0, 0, SRCCOPY)
       iretb  = endPaint( wnd%hWnd, ps )
       call LeaveCriticalSection(  loc(lpCriticalSection) )
     case (WM_LBUTTONDOWN)
       write(*, *) 'LBUTTON ', trim(wnd%title), wnd%hWnd
     case (WM_RBUTTONDOWN)
       iretb = DeleteObject( wnd%hDC )
       call PostQuitMessage( 0 )
     case default
       MainWndProc = DefWindowProc( hWnd, mesg, wParam, lParam )
    end select 
    return
  end function MainWndProc
!-------------------------------------------------------------------------------
  integer (HANDLE) function get_DosHndl()
    implicit none
    interface
      function GetConsoleWindow() ! non-existent in IFORTxx.MOD
      !DEC$ ATTRIBUTES DEFAULT, STDcall, DECORATE, ALIAS:'GetConsoleWindow' :: GetConsoleWindow
        use ifwinty
        integer (HANDLE) :: GetConsoleWindow
      end function
    end interface
    get_DosHndl = GetConsoleWindow()
    return
  end function get_DosHndl
!-------------------------------------------------------------------------------------
end module uho_win
!=================================================================================
module uhoplot
  use uho_win
  contains
!-------------------------------------------------------------------------------------
  subroutine gr_on(text, nx, ny)
    use IFMT ! multithread module
    implicit none
    character (LEN = *), intent(in), optional :: text
    integer           , intent(in), optional :: nx, ny
    integer (BOOL)    :: iretb
    integer (HANDLE)  :: hBmp
    integer (LPDWORD) :: id
    type (T_RECT)    :: rc
    if ( present(text) ) wnd%title   = TRIM(text)
    if ( present(nx)   ) wnd%nsize_x = nx
    if ( present(ny)   ) wnd%nsize_y = ny
    call InitializeCriticalSection( loc(lpCriticalSection))
    wnd%hThread = CreateThread(0, 0, Thread_Proc,  loc(wnd), CREATE_SUSPendED, id)
    wnd%id      = id
    iretb       = SetThreadPriority(wnd%hThread, THREAD_PRIORITY_BELOW_NORMAL)
    iretb       = ResumeThread(wnd%hThread)
    call sleep(100) ! wait for Window initialization 
    iretb = GetClientRect(wnd%hWnd, rc)
    hBmp  = CreateCompatibleBitmap(wnd%hDC, rc%right - rc%left, rc%bottom - rc%top)
    iretb = SelectObject(wnd%hDC, hBmp)
    iretb = DeleteObject(hBmp)
    iretb = PatBlt(wnd%hDC, 0, 0, rc%right - rc%left, rc%bottom - rc%top, WHITENESS)
    wnd%hPen = CreatePen(PS_SOLID, 1, 0)
    return
  end subroutine gr_on
!-------------------------------------------------------------------------------------
  subroutine gr_off(isec)
    use IFMT ! Module for multithread
    implicit none
    integer (DWORD), intent(in), optional :: isec 
    integer (BOOL)  :: iretb
    integer (DWORD) :: iwait
    if (present(isec)) THEN
      iwait = isec * 1000
    else
      iwait = INFINITE
    end if 
    call gr_show() 
    iretb = DeleteObject(wnd%hPen) 
    iretb = WaitForSingleObject(wnd%hThread, iwait)
    iretb = CloseHandle(wnd%hThread)
    iretb = PostMessage(wnd%hWnd, WM_DESTROY, NULL, NULL)
    wnd%hThread = NULL
    call DeleteCriticalSection( loc(lpCriticalSection) )
    return
  end subroutine gr_off
  !-------------------------------------------------------------------------------------
  integer (LONG) function Thread_Proc(lp_ThreadParameter)
  !DEC$ ATTRIBUTES STDcall, ALIAS:"_thread_proc" :: Thread_Proc
    implicit none
    integer (LPINT), intent(in) :: lp_ThreadParameter
    integer :: hInst
    type (t_wnd) :: wnd
    pointer (p_wnd, wnd) ! non-standard fortran
    p_wnd       = lp_ThreadParameter
    hInst       = GetModuleHandle(NULL) 
    Thread_Proc = WinMain(hInst, NULL, NULL, SW_SHOWNORMAL, wnd)
    return
  end function Thread_Proc
  !-------------------------------------------------------------------------------------
  subroutine gr_dot(ix, iy, icol)
    implicit none
    integer, intent(in) :: ix, iy
    integer, intent(in) :: icol
    integer (BOOL):: iretb
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb = SetPixel(wnd%hDC, ix, iy, icol)
    call LeaveCriticalSection( loc(lpCriticalSection) )
    return
  end subroutine gr_dot
  !-------------------------------------------------------------------------------------
  subroutine gr_move(ix, iy)
    implicit none
    integer, intent(in) :: ix, iy
    integer (BOOL):: iretb
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb = MoveToEx(wnd%hDC, ix, iy, NULL)
    call LeaveCriticalSection( loc(lpCriticalSection) )
    return
  end subroutine gr_move
  !-------------------------------------------------------------------------------------
  subroutine gr_line(ix, iy)
    implicit none
    integer, intent(in) :: ix, iy
    integer (BOOL):: iretb
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb = LineTo(wnd%hDC, ix, iy)
    call LeaveCriticalSection( loc(lpCriticalSection) )
    return
  end subroutine gr_line
  !-------------------------------------------------------------------------------------
  subroutine gr_show()
    implicit none
    integer (BOOL):: iretb
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb = InvalidateRect(wnd%hWnd, NULL, FALSE)
    call LeaveCriticalSection( loc(lpCriticalSection))
    return
  end subroutine gr_show
  !-------------------------------------------------------------------------------------
  subroutine gr_pencol(icol)
    implicit none
    integer, intent(in) :: icol
    integer :: iretb
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb    = DeleteObject(wnd%hPen) 
    wnd%hPen = CreatePen(PS_SOLID, 1, icol)
    iretb    = SelectObject(wnd%hDC, wnd%hPen)
    call LeaveCriticalSection( loc(lpCriticalSection))
    return
  end subroutine gr_pencol
  !-------------------------------------------------------------------------------------
  subroutine gr_text(ix, iy, txt, icol, ifontsize, ifontdirection)
    implicit none
    integer, intent(in) :: ix, iy
    character (LEN = *), intent(in) :: txt
    integer, intent(in), optional :: icol
    integer (BOOL)   :: iretb
    integer (HANDLE) :: hFont
    integer, optional, intent(in) :: ifontsize, ifontdirection
    integer :: kfontsize, kfontdirection

    if ( present(icol) ) iretb = SetTextColor(wnd%hDC, icol)
    if ( present(ifontsize) ) then 
      kfontsize = ifontsize
    else
      kfontsize = 10
    end if
    if ( present(ifontdirection) ) then
      kfontdirection = ifontdirection
    else
      kfontdirection = 0
    end if
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb = SetBkMode(wnd%hDC, TRANSPARENT)
    hFont = CreateFont( kfontsize , 10 , kfontdirection , 0 ,FW_DONTCARE , FALSE , FALSE , FALSE ,  &
		      	      ANSI_CHARSET , OUT_DEFAULT_PRECIS ,                   &
	  		          CLIP_DEFAULT_PRECIS , PROOF_QUALITY ,                 &
			          ior(FIXED_PITCH,  FF_ROMAN) , NULL       		)
    iretb = SelectObject(wnd%hdc , hFont)
    iretb = TextOut(wnd%hDC, ix, iy, txt, len_trim(txt))
    iretb = SelectObject(wnd%hdc , GetStockObject(SYSTEM_FONT))
    iretb = DeleteObject(hFont)
    call LeaveCriticalSection( loc(lpCriticalSection) )
    return
  end subroutine gr_text
  !-------------------------------------------------------------------------------------
  integer function irgb(ir, ig, ib)
    implicit none
    integer, intent(in) :: ir, ig, ib
    irgb = ir + (ig + (ib * 256)) * 256
    return
  end function irgb
!Relative coordinate  
  !-------------------------------------------------------------------------------------
  subroutine gr_axis(xmin, xmax, ymin, ymax)
    real, intent(in) :: xmin, xmax, ymin, ymax
    wnd%xmin = xmin
    wnd%xmax = xmax
    wnd%ymin = ymin
    wnd%ymax = ymax
    return
  end subroutine gr_axis
!----------------------------------------------------------------
  subroutine to_ixy(x, y, ix, iy)
    real   , intent(in ) :: x, y
    integer, intent(out) :: ix, iy 
    ix =  int( wnd%nsize_x * (x - wnd%xmin) / (wnd%xmax - wnd%xmin) )
    iy = -int( wnd%nsize_y * (y - wnd%ymin) / (wnd%ymax - wnd%ymin) ) + wnd%nsize_y
    return
  end subroutine to_ixy
!----------------------------------------------------------------
  subroutine gr_wdot(x, y, icol)
    real, intent(in) :: x, y
    integer, intent(in) :: icol
    integer :: ix, iy
    integer (BOOL):: iretb
    call to_ixy(x, y, ix, iy)
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb = SetPixel(wnd%hDC, ix, iy, icol)
    call LeaveCriticalSection( loc(lpCriticalSection) )
    return
  end subroutine gr_wdot
!----------------------------------------------------------------
  subroutine gr_wmove(x, y)
    real, intent(in) :: x, y
    integer (BOOL):: iretb
    integer :: ix, iy
    call to_ixy(x, y, ix, iy)
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb = MoveToEx(wnd%hDC, ix, iy, NULL)
    call LeaveCriticalSection( loc(lpCriticalSection) )
  return
  end subroutine gr_wmove
!----------------------------------------------------------------
  subroutine gr_wline(x, y)
    real, intent(in) :: x, y
    integer (BOOL):: iretb
    integer :: ix, iy
    call to_ixy(x, y, ix, iy)
    call EnterCriticalSection( loc(lpCriticalSection) )
    iretb = LineTo(wnd%hDC, ix, iy)
    call LeaveCriticalSection( loc(lpCriticalSection) )
    return
  end subroutine gr_wline
!----------------------------------------------------------------
  subroutine gr_wtext(x, y, text, icol, ifontsize, ifontdirection)
    real, intent(in) :: x, y
    character, intent(in) :: text*(*)
    integer, optional, intent(in) :: icol
    integer :: ix, iy
    integer (BOOL):: iretb
    integer, optional, intent(in) :: ifontsize, ifontdirection
    integer :: kfontsize, kfontdirection
    
    if ( present(icol) ) iretb = SetTextColor(wnd%hDC, icol)
    if ( present(ifontsize) ) then 
      kfontsize = ifontsize
    else
      kfontsize = 10
    end if
    if ( present(ifontdirection) ) then
      kfontdirection = ifontdirection
    else
      kfontdirection = 0
    end if
    call to_ixy(x, y, ix, iy)
    call EnterCriticalSection( loc(lpCriticalSection) )
    call gr_text(ix, iy, text, icol, kfontsize, kfontdirection)
    call LeaveCriticalSection( loc(lpCriticalSection) )
    return
  end subroutine gr_wtext
!----------------------------------------------------------------
  end module uhoplot
!=====================================================================================
module plotter
  use uhoplot
  real :: xpen = 0.0, ypen = 0.0
 contains
!----------------------------------------------------------------
  subroutine move(x, y)
  real, intent(in) :: x, y
    xpen = x
    ypen = y
    return
  end subroutine move
!----------------------------------------------------------------
  subroutine move_rel(dx, dy)
    real, intent(in) :: dx, dy
    xpen = xpen + dx
    ypen = ypen + dy
    return
  end subroutine move_rel
!----------------------------------------------------------------
  subroutine draw(x, y)
    real, intent(in) :: x, y
    call gr_wmove(xpen, ypen)
    call gr_wline(x   , y   )
    xpen = x
    ypen = y
    call gr_show()
    return
  end subroutine draw
!----------------------------------------------------------------
  subroutine draw_rel(dx, dy)
    real, intent(in) :: dx, dy
    call gr_wmove(xpen, ypen)
    call gr_wline(xpen + dx, ypen + dy)
    xpen = xpen + dx
    ypen = ypen + dy
    call gr_show()
    return
  end subroutine draw_rel
!----------------------------------------------------------------
end MODULE plotter