종합 실습 - 2차원 판 열 확산 시뮬레이터

이 장에서 배우는 것

열 확산 공식을 계산하는 것과 계산 프로그램을 완성하는 것은 서로 다른 일이다. 계산식이 맞아도 설정값을 잘못 읽으면 다른 문제를 풀게 된다. 저장한 온도 배열을 읽어도 시간과 계산 조건을 함께 복원하지 않으면 이어서 계산한 결과의 의미가 달라진다. 마지막으로 출력한 평균 온도가 경계까지 포함하는지 불분명하면 결과를 비교하기 어렵다.

이 장에서는 지금까지 만든 기능을 하나의 실행 흐름으로 묶는다. 작은 금속판에 초기 온도를 배치하고, 설정 파일에 지정한 횟수만큼 열 확산을 계산한다. 중간에는 체크포인트(checkpoint)를 저장하고 별도의 객체로 복원한다. 이후 계산을 마치면 계산 시각과 온도 통계를 일정한 서식으로 출력한다. 새로운 수치 기법을 추가하기보다 각 기능이 주고받는 상태와 책임을 명확히 하는 데 집중한다.

• 설정 읽기, 검증, 초기화, 계산, 저장, 복원, 요약 출력을 하나의 프로그램으로 연결한다.
• 격자와 시간 조건을 검증하여 계산을 시작하기 전에 잘못된 입력을 발견한다.
• 온도 배열과 계산 조건을 함께 저장하여 재시작 상태의 일관성을 확인한다.
• 직렬 실행과 OpenMP 실행에서 같은 격자 갱신 절차를 사용한다.
• 작은 기준 사례의 예상값으로 완성 프로그램의 동작을 확인한다.

문제 상황

얇은 금속판의 가운데를 잠시 가열한 뒤 열이 주변으로 퍼지는 과정을 계산한다고 하자. 판 가장자리는 냉각 장치에 붙어 있어 항상 온도 0으로 유지된다. 판 내부에는 열원이 더 이상 없으며, 열 확산 계수는 위치와 시간에 따라 변하지 않는다. 여기서 온도 0은 기준 온도와의 차이가 0이라는 뜻으로 해석할 수 있다.

개발 중에는 계산을 몇 번만 수행하지만 실제 작업에서는 계산 조건을 바꾸며 반복 실행하게 된다. 계산을 중단했다가 이어서 실행할 수도 있다. 이때 배열만 저장하는 것으로는 충분하지 않다. 같은 배열이라도 격자 간격이나 시간 간격이 다르면 다음 단계의 온도가 달라진다. 계산 결과에는 온도뿐 아니라 그 온도를 만든 조건이 따라다녀야 한다.

이번 예제는 재시작 과정을 짧은 실행 안에서 확인하도록 구성한다. 기본 설정으로 두 단계를 계산한 뒤 파일을 저장한다. 이어서 새 판 객체가 파일을 읽고 나머지 두 단계를 계산한다. 운영 프로그램에서는 저장 직후 종료하고 다른 실행에서 복원할 수 있지만, 여기서는 저장 전후의 연결을 눈으로 확인할 수 있게 한 실행 안에 두 경로를 배치한다.

격자는 경계를 포함하여 5행 5열이다. 내부 3행 3열의 가운데에만 초기 온도 100을 넣는다. 작은 격자를 선택한 이유는 결과를 직접 계산해 검증하기 위해서다. 화면에 보이는 출력은 작아도 설정 검증, 객체 상태 관리, 파일 오류 확인은 더 큰 계산에서도 필요한 구조를 따른다.

실행 흐름과 책임을 나누기

프로그램은 설정 타입과 판 타입을 가진 모듈 하나, 그리고 실행 순서를 정하는 주 프로그램으로 구성한다. 설정 타입은 계산 조건을 보관한다. 판 타입은 온도 배열, 다음 온도 배열, 완료한 단계 수, 현재 계산 시각을 보관한다. 주 프로그램은 언제 저장하고 언제 복원할지 결정한다. 판 객체가 파일 이름이나 전체 실행 횟수를 임의로 정하지 않도록 책임을 나눈다.

완성 프로그램에서 각 부분이 맡는 책임

부분	보관하는 정보	수행하는 일
설정 타입	격자 크기, 계산 횟수, 물리 계수	계산 조건을 하나의 값으로 전달한다.
판 타입	설정, 두 배열, 단계 수, 시각	초기화, 한 단계 계산, 저장, 복원을 수행한다.
주 프로그램	실행 중 사용할 객체 두 개	저장 시점과 계산 종료 시점을 정한다.
검사 프로시저	별도 상태 없음	조건이 맞지 않으면 원인을 표시하고 종료한다.

설정 파일은 이름 목록(namelist)으로 읽는다. 기본서를 마친 독자에게 필요한 핵심은 문법 자체보다 읽은 값의 검증이다. 이름 목록은 항목 일부만 적을 수 있으므로, 읽기 전에 모든 항목에 기본값을 넣는다. 그러면 파일에 생략된 항목은 기본값을 유지한다. 처음 실행할 때 파일이 없으면 프로그램이 기본 파일을 만든다. 이미 있는 파일은 덮어쓰지 않으므로 사용자가 다음 실행의 조건을 수정할 수 있다.

설정 검증은 초기화 전에 끝낸다. 격자가 너무 작으면 내부점이 없고, 시간 간격이 음수면 계산 방향이 의도와 달라진다. 실수에 무한대나 비수가 들어오면 양수 비교만으로 충분히 걸러지지 않을 수 있다. 따라서 유한한 값인지 먼저 검사한 뒤 양수 조건과 안정성 조건을 확인한다. 여기서 허용한 격자 크기와 계산 횟수의 상한은 예제의 작업 범위를 정하는 제한이며 Fortran 언어의 제한은 아니다.

주 프로그램의 두 객체는 역할이 다르다. 첫 번째 객체는 초기 상태에서 저장 시점까지 계산한다. 두 번째 객체는 저장 파일에서 출발한다. 두 객체를 사용하면 복원이 실제로 수행됐는지 확인하기 쉽다. 기존 객체의 배열이 우연히 남아 있어 파일을 제대로 읽지 않았는데도 이어서 계산되는 상황을 피할 수 있다.

한 단계 계산과 경계의 계약

앞에서 사용한 명시적 유한차분식을 그대로 적용한다. x 방향 계수를 rx = alpha × dt / dx², y 방향 계수를 ry = alpha × dt / dy²로 둔다. 내부점의 다음 온도는 현재 온도에 두 방향의 이웃 온도 차이를 더해 구한다. 이 방법에서는 rx + ry가 0.5 이하가 되도록 시간 간격을 선택한다. 프로그램은 이 조건을 확인한 뒤 계산을 시작한다.

한 단계 계산의 중요한 계약은 현재 온도를 읽는 동안 그 배열을 바꾸지 않는다는 것이다. 모든 내부점의 새 온도를 보조 배열에 쓴 다음, 갱신이 끝나면 보조 배열 전체를 현재 배열로 복사한다. 이 구조를 유지하면 내부점 계산을 여러 스레드에 나누어도 각 반복이 서로 다른 위치에 기록한다. 반복 사이에 합계나 최댓값을 누적하지도 않는다.

경계는 매 단계 보조 배열을 0으로 채우는 동작으로 유지한다. 내부점만 갱신하기 때문에 가장자리 네 변은 계속 0이다. 복원할 때도 파일 속 경계가 0인지 검사한다. 경계를 고정한다는 가정은 계산식 바깥의 장식이 아니라 이 문제를 정의하는 조건이다.

기본 설정에서는 alpha가 0.1이고 dt, dx, dy가 모두 1이므로 rx와 ry는 각각 0.1이다. 초기 가운데 온도 100은 한 단계 뒤 60이 되고, 바로 인접한 네 점은 각각 10이 된다. 경계와 접한 내부점에 열이 도달하면 열이 경계를 통해 빠져나간다. 따라서 내부 온도 합계가 계속 100으로 유지되는 것을 기대하면 안 된다.

예제는 바깥쪽 j 반복에 OpenMP 지시문을 붙인다. -fopenmp 옵션이 없으면 지시문은 주석으로 취급되어 같은 반복문이 직렬로 실행된다. -fopenmp를 사용하면 격자 행을 나누어 계산한다. 각 점에서 수행하는 산술 순서는 바꾸지 않고, 요약 통계는 병렬 구간 바깥에서 계산한다. 스레드 수나 실행 시간은 출력하지 않아 기준 출력이 실행 환경에 따라 달라지지 않게 한다.

do concurrent도 반복 사이의 독립성을 표현하는 방법이다. 이번 완성 코드는 스레드 실행을 명시적으로 선택할 수 있도록 OpenMP 반복문을 채택한다. 같은 반복에 두 표현을 무조건 겹쳐 넣기보다, 현재 배열을 읽고 보조 배열만 쓰는 구조가 먼저 성립하는지 확인하는 것이 중요하다.

재시작 파일에 담아야 하는 상태

체크포인트에는 파일 식별 문자열, 격자 크기, 완료 단계, 계산 시각, 열 확산 계수, 시간 간격, 두 격자 간격, 온도 배열을 순서대로 기록한다. 식별 문자열 PLATE001의 끝 숫자는 이 예제의 파일 형식을 구분한다. 기록 순서를 바꾸는 경우에는 식별 문자열도 바꾸어 이전 파일을 새 형식으로 잘못 읽지 않게 해야 한다.

저장하는 상태와 다시 만드는 상태의 구분

정보	저장 여부	복원할 때의 처리
현재 온도 배열	저장한다.	유한한 값과 경계 조건을 검사한다.
단계 수와 계산 시각	저장한다.	단계 범위와 시각의 일관성을 검사한다.
격자와 계산 계수	저장한다.	현재 설정과 일치하는지 확인한다.
다음 온도 배열	저장하지 않는다.	배열을 새로 할당하고 다음 단계에서 덮어쓴다.

보조 배열은 다음 계산을 위한 작업 공간이므로 저장하지 않는다. 전체 종료 단계 수도 파일에 묶지 않는다. 저장 당시보다 더 오래 계산하고 싶을 수 있기 때문이다. 대신 복원된 완료 단계가 현재 설정의 종료 단계보다 크지 않은지 확인한다. 저장 시점을 정하는 설정 역시 온도 상태의 해석에는 필요하지 않다.

복원은 임시 판 객체에 먼저 수행한다. 파일의 머리 부분을 검사하고, 크기에 맞는 배열을 준비한 다음 온도를 읽는다. 온도와 시각 검사를 모두 통과해야 대상 객체에 대입한다. 파생 타입에 포함된 할당 배열은 이 대입에서 값이 복사된다. 임시 객체가 사라져도 대상 객체가 같은 메모리를 가리키는 구조가 아니다.

이 파일은 같은 프로그램과 호환되는 실행 환경에서 재시작하기 위한 스트림 이진 파일이다. 기본 정수의 크기, 실수 표현, 바이트 순서까지 지정한 교환 형식은 아니다. 또한 파일 식별과 값 검사는 손상 검사의 일부일 뿐이다. 내부 온도 한 개가 다른 유한한 값으로 바뀐 손상까지 모두 찾아내지는 못한다. 여러 환경 사이의 교환이나 강한 손상 검사가 필요하면 별도의 형식 규약과 검증 정보를 설계해야 한다.

완성 코드

다음 내용을 main.f90으로 저장한다. 모듈을 같은 파일 앞부분에 두었으므로 별도의 소스 파일은 필요하지 않다. settings.dat가 없는 디렉터리에서 실행하면 프로그램이 먼저 기본 설정 파일을 만든다. 표준 입력은 사용하지 않는다.

main.f90

module heat_project
  use iso_fortran_env, only: real64
  use ieee_arithmetic, only: ieee_is_finite
  implicit none
  private

  public :: settings_t, plate_t, read_settings, require

  type :: settings_t
    integer :: nx = 5, ny = 5
    integer :: steps = 4, checkpoint_step = 2
    real(real64) :: alpha = 0.1_real64
    real(real64) :: dt = 1.0_real64
    real(real64) :: dx = 1.0_real64, dy = 1.0_real64
  end type settings_t

  type :: plate_t
    type(settings_t) :: cfg
    real(real64), allocatable :: temperature(:, :)
    real(real64), allocatable :: work(:, :)
    integer :: completed = 0
    real(real64) :: time = 0.0_real64
  contains
    procedure :: initialize
    procedure :: advance
    procedure :: save
    procedure :: restore
    procedure :: interior_mean
  end type plate_t

contains

  subroutine require(condition, message)
    logical, intent(in) :: condition
    character(*), intent(in) :: message

    if (.not. condition) error stop message
  end subroutine require

  subroutine read_settings(filename, cfg)
    character(*), intent(in) :: filename
    type(settings_t), intent(out) :: cfg
    integer :: nx, ny, steps, checkpoint_step
    integer :: unit, status
    real(real64) :: alpha, dt, dx, dy, rx, ry
    logical :: exists
    namelist /run/ nx, ny, steps, checkpoint_step, alpha, dt, dx, dy

    cfg = settings_t()
    nx = cfg%nx
    ny = cfg%ny
    steps = cfg%steps
    checkpoint_step = cfg%checkpoint_step
    alpha = cfg%alpha
    dt = cfg%dt
    dx = cfg%dx
    dy = cfg%dy

    inquire(file=filename, exist=exists, iostat=status)
    call require(status == 0, 'cannot inspect settings file')
    if (.not. exists) then
      open(newunit=unit, file=filename, status='new', &
           action='write', iostat=status)
      call require(status == 0, 'cannot create settings file')
      write(unit, '(a)', iostat=status) &
        '&run', &
        ' nx=5, ny=5, steps=4, checkpoint_step=2,', &
        ' alpha=0.1, dt=1.0, dx=1.0, dy=1.0', &
        '/'
      call require(status == 0, 'cannot write settings file')
      close(unit, iostat=status)
      call require(status == 0, 'cannot close settings file')
    end if

    open(newunit=unit, file=filename, status='old', &
         action='read', iostat=status)
    call require(status == 0, 'cannot open settings file')
    read(unit, nml=run, iostat=status)
    call require(status == 0, 'invalid settings file')
    close(unit, iostat=status)
    call require(status == 0, 'cannot close settings file')

    call require(nx >= 3 .and. nx <= 501, 'invalid nx')
    call require(ny >= 3 .and. ny <= 501, 'invalid ny')
    call require(steps >= 1 .and. steps <= 1000000, &
                 'invalid steps')
    call require(checkpoint_step >= 1 .and. &
                 checkpoint_step <= steps, 'invalid checkpoint step')
    call require(all(ieee_is_finite([alpha, dt, dx, dy])), &
                 'nonfinite settings')
    call require(min(alpha, dt, dx, dy) > 0.0_real64, &
                 'settings must be positive')

    rx = alpha * dt / dx**2
    ry = alpha * dt / dy**2
    call require(all(ieee_is_finite([rx, ry])), &
                 'nonfinite diffusion weights')
    call require(rx > 0.0_real64 .and. ry > 0.0_real64, &
                 'diffusion weights must be positive')
    call require(rx + ry <= 0.5_real64, 'unstable time step')
    call require(ieee_is_finite(real(steps, real64) * dt), &
                 'nonfinite final time')

    cfg = settings_t(nx, ny, steps, checkpoint_step, alpha, dt, dx, dy)
  end subroutine read_settings

  subroutine initialize(this, cfg)
    class(plate_t), intent(inout) :: this
    type(settings_t), intent(in) :: cfg
    integer :: status

    this%cfg = cfg
    if (allocated(this%temperature)) deallocate(this%temperature)
    if (allocated(this%work)) deallocate(this%work)
    allocate(this%temperature(cfg%nx, cfg%ny), &
             this%work(cfg%nx, cfg%ny), stat=status)
    call require(status == 0, 'cannot allocate plate')

    this%temperature = 0.0_real64
    this%work = 0.0_real64
    this%temperature((cfg%nx + 1) / 2, (cfg%ny + 1) / 2) = 100.0_real64
    this%completed = 0
    this%time = 0.0_real64
  end subroutine initialize

  subroutine advance(this)
    class(plate_t), intent(inout) :: this
    integer :: i, j
    real(real64) :: rx, ry

    call require(this%completed < this%cfg%steps, &
                 'requested step exceeds final step')
    rx = this%cfg%alpha * this%cfg%dt / this%cfg%dx**2
    ry = this%cfg%alpha * this%cfg%dt / this%cfg%dy**2
    this%work = 0.0_real64

    !$omp parallel do default(none) shared(this, rx, ry) private(i, j)
    do j = 2, this%cfg%ny - 1
      do i = 2, this%cfg%nx - 1
        this%work(i, j) = this%temperature(i, j) &
          + rx * (this%temperature(i - 1, j) &
          - 2.0_real64 * this%temperature(i, j) &
          + this%temperature(i + 1, j)) &
          + ry * (this%temperature(i, j - 1) &
          - 2.0_real64 * this%temperature(i, j) &
          + this%temperature(i, j + 1))
      end do
    end do
    !$omp end parallel do

    this%temperature = this%work
    this%completed = this%completed + 1
    this%time = real(this%completed, real64) * this%cfg%dt
  end subroutine advance

  subroutine save(this, filename)
    class(plate_t), intent(in) :: this
    character(*), intent(in) :: filename
    integer :: unit, status

    open(newunit=unit, file=filename, status='replace', &
         access='stream', form='unformatted', action='write', &
         iostat=status)
    call require(status == 0, 'cannot create checkpoint')
    write(unit, iostat=status) 'PLATE001', &
      this%cfg%nx, this%cfg%ny, this%completed, this%time, &
      this%cfg%alpha, this%cfg%dt, this%cfg%dx, this%cfg%dy, &
      this%temperature
    call require(status == 0, 'cannot write checkpoint')
    close(unit, iostat=status)
    call require(status == 0, 'cannot close checkpoint')
  end subroutine save

  subroutine restore(this, filename, cfg)
    class(plate_t), intent(inout) :: this
    character(*), intent(in) :: filename
    type(settings_t), intent(in) :: cfg
    type(plate_t) :: candidate
    character(8) :: signature
    integer :: unit, status, nx, ny, completed
    real(real64) :: saved_time, alpha, dt, dx, dy, expected_time

    open(newunit=unit, file=filename, status='old', &
         access='stream', form='unformatted', action='read', &
         iostat=status)
    call require(status == 0, 'cannot open checkpoint')
    read(unit, iostat=status) signature, nx, ny, completed, &
                             saved_time, alpha, dt, dx, dy
    call require(status == 0, 'cannot read checkpoint header')
    call require(signature == 'PLATE001', 'unknown checkpoint format')
    call require(nx == cfg%nx .and. ny == cfg%ny, &
                 'checkpoint grid mismatch')
    call require(completed >= 0 .and. completed <= cfg%steps, &
                 'checkpoint step out of range')
    call require(all(ieee_is_finite([saved_time, alpha, dt, dx, dy])), &
                 'nonfinite checkpoint header')
    call require(maxval(abs([alpha - cfg%alpha, dt - cfg%dt, &
                 dx - cfg%dx, dy - cfg%dy])) <= 0.0_real64, &
                 'checkpoint settings mismatch')

    expected_time = real(completed, real64) * cfg%dt
    call require(abs(saved_time - expected_time) <= &
                 8.0_real64 * epsilon(1.0_real64) &
                 * max(1.0_real64, expected_time), &
                 'checkpoint time mismatch')

    call candidate%initialize(cfg)
    read(unit, iostat=status) candidate%temperature
    call require(status == 0, 'cannot read checkpoint temperatures')
    close(unit, iostat=status)
    call require(status == 0, 'cannot close checkpoint')
    call require(all(ieee_is_finite(candidate%temperature)), &
                 'nonfinite checkpoint temperatures')
    call require(maxval(abs(candidate%temperature(1, :))) <= &
                 0.0_real64 .and. &
                 maxval(abs(candidate%temperature(nx, :))) <= &
                 0.0_real64 .and. &
                 maxval(abs(candidate%temperature(:, 1))) <= &
                 0.0_real64 .and. &
                 maxval(abs(candidate%temperature(:, ny))) <= &
                 0.0_real64, 'checkpoint boundary mismatch')

    candidate%completed = completed
    candidate%time = expected_time
    select type (this)
    type is (plate_t)
      this = candidate
    class default
      error stop 'unsupported extended plate type'
    end select
  end subroutine restore

  function interior_mean(this) result(value)
    class(plate_t), intent(in) :: this
    real(real64) :: value
    integer :: nx, ny

    nx = this%cfg%nx
    ny = this%cfg%ny
    value = sum(this%temperature(2:nx - 1, 2:ny - 1)) &
            / real((nx - 2) * (ny - 2), real64)
  end function interior_mean

end module heat_project

program main
  use heat_project, only: settings_t, plate_t, read_settings
  implicit none
  type(settings_t) :: cfg
  type(plate_t) :: original, resumed
  integer :: cx, cy

  call read_settings('settings.dat', cfg)
  call original%initialize(cfg)
  do while (original%completed < cfg%checkpoint_step)
    call original%advance()
  end do

  call original%save('plate.chk')
  call resumed%restore('plate.chk', cfg)
  do while (resumed%completed < cfg%steps)
    call resumed%advance()
  end do

  cx = (cfg%nx + 1) / 2
  cy = (cfg%ny + 1) / 2
  write(*, '(a,i0,a,i0)') 'grid: ', cfg%nx, ' x ', cfg%ny
  write(*, '(a,i0)') 'checkpoint step: ', original%completed
  write(*, '(a,i0)') 'completed steps: ', resumed%completed
  write(*, '(a,f10.4)') 'time: ', resumed%time
  write(*, '(a,f10.4)') 'center: ', resumed%temperature(cx, cy)
  write(*, '(a,f10.4)') 'interior mean: ', resumed%interior_mean()
  write(*, '(a,f10.4)') 'maximum: ', maxval(resumed%temperature)
end program main

줄별 해설

모듈 첫 부분의 private는 내부 구현을 기본적으로 감춘다. public 목록에 설정 타입, 판 타입, 설정 읽기와 검사 프로시저만 공개한다. real64를 가져와 배열, 계산 계수, 저장 파일의 실수 종류를 통일한다. 파일에서 읽은 실수도 같은 종류로 선언해야 기록한 바이트를 의도한 값으로 해석할 수 있다.

settings_t의 성분 기본값은 설정을 읽기 전의 출발점이다. cfg = settings_t()로 기본 상태를 만든 뒤 이름 목록에 들어갈 지역 변수로 옮긴다. 이름 목록 입력에는 판 객체 전체를 노출하지 않는다. 파일에서 바꿀 수 있는 항목을 제한하고, 검증을 통과한 값만 마지막 구조 생성자에서 cfg에 넣는다. 구조 생성자의 위치 인수는 타입 성분 선언 순서와 맞아야 한다.

파일을 만드는 open에는 status='new'를 사용한다. 파일이 없다는 확인 뒤에도 다른 작업이 같은 이름의 파일을 만들었다면 덮어쓰지 않고 실패하게 된다. 반면 체크포인트 저장은 이번 실행의 저장 상태로 교체하려는 의도이므로 status='replace'를 사용한다. 이 예제는 하나의 체크포인트만 유지하며, 교체 도중의 장애에 대비하는 임시 파일 방식까지 구현하지는 않는다.

require는 검사 조건을 한곳에서 표현하게 한다. 파일을 열 수 없을 때, 입력 문법이 틀렸을 때, 저장 파일의 조건이 다를 때 서로 다른 메시지가 나온다. iostat를 모든 주요 입출력에 붙이고 close의 결과까지 확인한다. 기록 요청이 성공했더라도 파일을 닫는 단계에서 오류가 발생할 수 있기 때문이다. 정상 실행 결과와 오류 종료 때 런타임이 덧붙이는 진단 출력은 구분해서 읽는다.

initialize는 기존 배열이 있으면 해제한 뒤 두 배열을 함께 할당한다. 현재 온도와 작업 배열을 모두 0으로 채우고 가운데 한 점만 100으로 바꾼다. 짝수 격자에서는 가운데가 한 점으로 정해지지 않으므로 이 코드의 정수 나눗셈은 중앙의 두 인덱스 중 작은 쪽을 택한다. 기본 설정은 홀수 격자이므로 가운데 좌표가 명확하다.

advance의 시작 검사는 종료 단계 이후의 추가 호출을 막는다. 바깥 반복은 j, 안쪽 반복은 i로 두어 첫 번째 첨자가 연속으로 변하게 한다. OpenMP 구간에서 this는 공유되지만 읽는 temperature와 쓰는 work의 역할이 분리되어 있다. i와 j는 스레드별 반복 변수다. 구간 끝에서 모든 행의 계산이 끝난 후 배열 대입을 수행하므로 일부 행만 갱신된 상태를 다음 단계가 읽지 않는다.

시각은 이전 시각에 dt를 반복해서 더하지 않고 완료 단계에 dt를 곱해 구한다. 완료 단계가 재시작의 기준이 되며, 시각을 만드는 방법도 저장 전후에 동일해진다. 이것이 모든 시간 오차를 없애는 것은 아니다. 실수 곱셈의 반올림은 남지만, 단계 수와 시각을 연결하는 규칙을 하나로 유지할 수 있다.

save의 출력 목록과 restore의 입력 목록은 같은 순서를 따른다. 스트림 파일에는 이 목록의 항목 이름이 기록되지 않는다. 읽는 쪽이 정확한 종류와 순서를 알고 있어야 한다. 복원에서는 먼저 작은 머리 정보를 읽고 격자 일치를 확인한 다음 배열을 준비한다. 파일에서 받은 크기를 검사 없이 할당 크기로 사용하지 않는다.

설정 계수 비교는 저장 당시와 현재 값의 정확한 일치를 요구한다. 아주 작은 차이도 다른 계산 조건으로 취급한다. 반면 시각 검사는 반올림 규모의 허용 범위를 둔다. 시각은 단계 수에서 다시 계산할 수 있는 파생 값이므로 검사를 통과한 뒤 expected_time으로 정규화한다. 조건 계수와 파생 값에 같은 비교 정책을 기계적으로 적용하지 않는다는 점을 기억한다.

마지막 select type은 다형적인 통과 객체를 구체적인 plate_t로 좁혀 전체 객체 대입을 수행한다. 여기서는 확장 타입의 추가 상태를 어떻게 복원할지 정의하지 않았으므로 확장 객체는 거부한다. 확장을 허용하려면 추가 성분의 저장 규칙과 복원 계약도 함께 설계해야 한다. 단순히 기본 타입 부분만 복사하고 복원이 끝났다고 판단하면 확장 성분은 이전 상태로 남을 수 있다.

interior_mean은 경계가 빠진 배열 구간만 합산한다. 격자 크기 검증 덕분에 분모가 0이 되지 않는다. 주 프로그램은 저장 전 객체의 완료 단계를 출력하고, 최종 시각과 온도는 복원 후 객체에서 출력한다. 두 경로의 정보를 구분하면 체크포인트를 어느 시점에 만들고 어디까지 이어 계산했는지가 드러난다.

실행 결과

settings.dat가 없는 새 작업 디렉터리에서 다음 명령을 실행한다. 정상 컴파일에서는 경고가 나오지 않는 코드다. 실행하면 기본 설정 파일과 plate.chk가 생성된다.

gfortran -std=f2018 -Wall main.f90 -o plate
./plate

grid: 5 x 5
checkpoint step: 2
completed steps: 4
time:     4.0000
center:    21.9200
interior mean:     9.1200
maximum:    21.9200

OpenMP 실행은 다음과 같이 선택한다. 기본 설정에서는 출력의 마지막 네 자리까지 위 결과와 같다. 이 격자는 매우 작아 병렬화의 속도 이득을 평가하는 용도로는 적합하지 않다.

gfortran -std=f2018 -Wall -fopenmp main.f90 -o plate
OMP_NUM_THREADS=2 ./plate

프로그램이 만든 설정 파일의 내용은 다음과 같다. 이미 파일이 있으면 그 내용을 읽으므로, 아래 예상 결과를 다시 확인하려면 같은 설정을 사용해야 한다.

settings.dat

&run
 nx=5, ny=5, steps=4, checkpoint_step=2,
 alpha=0.1, dt=1.0, dx=1.0, dy=1.0
/

최종 내부 온도는 아래 배열과 같다. 행은 위에서 아래로 y 방향 위치를 표시한다. 이 사례는 두 방향 조건이 같고 초기 상태가 대칭이므로 x와 y를 바꾸어도 같은 배열이 된다. 경계를 포함한 전체 평균과 아래 내부 평균은 서로 다르다.

    4.4800   10.5600    4.4800
   10.5600   21.9200   10.5600
    4.4800   10.5600    4.4800

가운데 값 21.92, 인접한 네 점의 합 42.24, 모서리 네 점의 합 17.92를 더하면 내부 합계는 82.08이다. 이를 내부점 9개로 나누면 평균 9.12가 된다. 이러한 작은 기준 사례는 파일 복원, 반복 횟수, 경계 처리의 연결이 맞는지 확인하는 데 유용하다. 일반 조건에서의 정확도 검증을 이 한 사례로 대신할 수는 없다.

실무에서 자주 틀리는 것

현재 배열에 곧바로 새 온도를 쓰기

다음 코드는 앞서 갱신한 이웃을 뒤쪽 점이 읽게 한다. 반복 순서가 수치 방법의 일부가 되어, 의도한 명시적 갱신과 다른 계산이 된다. 병렬 실행에서도 읽기와 쓰기가 충돌할 수 있다.

! 잘못된 코드
this%temperature(i, j) = this%temperature(i, j) &
  + rx * (this%temperature(i - 1, j) &
  - 2.0_real64 * this%temperature(i, j) &
  + this%temperature(i + 1, j))

고친 코드는 현재 배열만 읽고 작업 배열에 쓴다. y 방향 항도 같은 방식으로 더하고, 전체 반복이 끝난 뒤 현재 배열을 교체한다.

! 고친 코드의 x 방향 부분
this%work(i, j) = this%temperature(i, j) &
  + rx * (this%temperature(i - 1, j) &
  - 2.0_real64 * this%temperature(i, j) &
  + this%temperature(i + 1, j))
! 전체 격자 계산이 끝난 뒤 수행한다.
this%temperature = this%work

복원 후 전체 횟수를 다시 계산하기

완료 단계가 2이고 종료 단계가 4인데 네 번을 더 계산하면 총 여섯 단계가 된다. 배열 복원이 성공해도 반복 범위를 잘못 잡으면 요청한 시각을 넘는다. 완성 코드의 advance는 이런 추가 호출을 검사하여 종료한다.

! 잘못된 코드
do k = 1, cfg%steps
  call resumed%advance()
end do

남은 계산은 저장된 완료 단계와 목표 단계의 차이로 정한다. 완료 단계를 객체가 관리하게 하면 주 프로그램에 별도의 단계 계수를 중복해서 둘 필요가 없다.

! 고친 코드
do while (resumed%completed < cfg%steps)
  call resumed%advance()
end do

온도 배열만 읽고 계산 조건을 바꾸기

다음 입력은 파일에 배열만 있다는 가정 아래에서는 읽을 수 있지만, 그 배열이 어떤 간격으로 계산됐는지 확인하지 않는다. 완성 코드의 파일에는 머리 정보가 있으므로 이런 방식으로 읽으면 기록 위치부터 맞지 않는다.

! 잘못된 코드
read(unit) resumed%temperature
resumed%cfg%dt = 0.5_real64

같은 문제의 연속 계산이 목적이라면 머리 정보와 현재 설정의 일치를 확인한다. 시간 간격 변경을 허용하는 별도 기능을 만들려면 완료 단계와 시각의 관계부터 다시 정의해야 한다.

! 고친 코드
call resumed%restore('plate.chk', cfg)
! restore 안에서 격자와 alpha, dt, dx, dy를 검사한다.

평균의 대상 영역을 설명하지 않기

전체 배열의 평균을 내부 평균이라고 표시하면 냉각 경계의 0이 분모에 포함된다. 기본 사례에서 전체 평균은 82.08을 25로 나눈 3.2832다. 계산값 자체보다 통계 이름과 대상 영역의 불일치가 문제다.

! 잘못된 코드: 내부 평균이라는 이름과 맞지 않는다.
value = sum(this%temperature) &
        / real(size(this%temperature), real64)

내부 평균을 원하면 합계와 개수 모두 같은 내부 영역을 기준으로 계산한다. 전체 평균이 필요하다면 기존 계산을 유지하고 출력 이름을 전체 평균으로 바꾼다.

! 고친 코드
value = sum(this%temperature(2:nx - 1, 2:ny - 1)) &
        / real((nx - 2) * (ny - 2), real64)

한눈에 보기

실행 단계마다 유지해야 하는 조건

단계	입력	확인할 조건	결과
설정 읽기	settings.dat	유한한 양수와 안정성 조건	검증된 설정 객체
초기화	설정 객체	내부점이 있는 격자	초기 온도와 작업 배열
한 단계 계산	현재 온도	현재 배열은 읽기 전용으로 사용	다음 온도와 증가한 단계
저장	판 객체	머리 정보와 배열의 기록 순서	plate.chk
복원	체크포인트와 현재 설정	조건 일치, 유한한 온도, 고정 경계	이어서 계산할 판 객체
요약	최종 판 객체	영역 정의와 고정 소수점 서식	비교 가능한 출력

완성 프로그램의 중심은 기능의 개수가 아니라 상태의 일관성이다. 설정은 계산 전에 검증하고, 한 단계 안에서는 읽는 배열과 쓰는 배열을 나누며, 파일을 읽은 상태는 검사를 통과한 뒤 채택한다. 출력에서도 단계, 시각, 통계 영역을 함께 명시한다. 이 연결이 유지되어야 계산 결과를 다시 실행하고 비교할 수 있다.

연습 문제

1. 기본 설정에서 steps를 2로 바꾸고 checkpoint_step은 2로 유지하라. 복원 후 추가 계산이 수행되는지 설명하고, 가운데 온도와 내부 평균의 예상값을 구하라.
2. 기본 설정의 dt를 3.0으로 바꾸어라. 프로그램이 계산을 시작할 수 있는지 rx와 ry를 이용해 판단하라. 다른 조건을 유지할 때 허용되는 dt의 상한도 구하라.
3. 저장과 복원을 거친 결과를 처음부터 끝까지 연속 계산한 결과와 비교하는 확인 코드를 추가하라. 별도 판 객체를 사용하고, 두 온도 배열의 최대 절대 차이를 고정 소수점 서식으로 출력하라.
4. 프로그램이 plate.chk를 먼저 찾고, 파일이 있으면 복원하여 계산하며 없으면 새 계산을 시작하도록 실행 경로를 바꾸어라. 기본 설정에서 두 번째 실행이 어느 단계부터 계산하는지 설명하라.

정답과 해설

1. 복원한 완료 단계가 종료 단계와 같은 2이므로 뒤쪽 do while은 한 번도 실행되지 않는다. 두 단계 뒤 가운데는 40.0000, 인접한 네 점은 각각 12.0000, 내부 모서리 네 점은 각각 2.0000이다. 내부 합계는 96이며 평균은 10.6667로 출력된다. 저장 시점과 종료 시점이 같은 경우도 별도 예외 분기 없이 처리된다.

2. dx와 dy가 1이고 alpha가 0.1이므로 두 계수는 각각 0.3이다. 합 0.6이 0.5보다 커서 unstable time step 메시지로 종료한다. 합은 0.2 × dt이므로 이 설정의 dt 상한은 2.5다. 컴파일 성공 여부와 수치 계산 조건의 적절함은 별개다.

3. 주 프로그램에서 모든 선언이 끝나기 전에 다음 선언을 추가한다. real64를 사용하려면 주 프로그램의 use 문 위치에 iso_fortran_env 가져오기도 추가한다.

   use iso_fortran_env, only: real64
   ! 아래 선언은 기존 선언들과 함께 둔다.
   type(plate_t) :: reference
   real(real64) :: difference

   기존 최종 계산이 끝난 뒤 다음 코드를 실행한다. 기본 사례의 출력은 difference: 뒤에 공백 다섯 개와 0.0000이 온다. 표시 자릿수가 제한되어 있으므로 출력이 0이라고 해서 임의의 두 결과가 비트 단위로 같다는 뜻은 아니다.

   call reference%initialize(cfg)
   do while (reference%completed < cfg%steps)
     call reference%advance()
   end do
   difference = maxval(abs(reference%temperature - resumed%temperature))
   write(*, '(a,f10.4)') 'difference: ', difference

   이 비교는 같은 계산식에서 재시작 경로가 상태를 보존하는지 확인한다. 연속 계산과 재시작 계산이 같은 오류를 공유할 수 있으므로 수치 방법 자체의 검증에는 별도의 기준값도 필요하다.

4. 주 프로그램에 논리 변수 exists와 정수 변수 status를 선언하고 require를 가져온다. 기존 초기화부터 최종 계산까지의 부분을 다음 흐름으로 바꾼다. 요약 출력에서는 original 대신 실제 계산 객체 resumed의 상태를 사용하고, 체크포인트 단계의 별도 출력은 삭제하거나 저장 정책에 맞게 다시 정의한다.

   inquire(file='plate.chk', exist=exists, iostat=status)
   call require(status == 0, 'cannot inspect checkpoint')
   if (exists) then
     call resumed%restore('plate.chk', cfg)
   else
     call resumed%initialize(cfg)
   end if
   
   do while (resumed%completed < cfg%steps)
     call resumed%advance()
     if (resumed%completed == cfg%checkpoint_step) then
       call resumed%save('plate.chk')
     end if
   end do

   처음 실행은 단계 0에서 시작하고 단계 2의 파일을 남긴 뒤 단계 4까지 계산한다. 두 번째 실행은 저장된 단계 2부터 출발하여 단계 3과 4를 계산한다. 여기서는 최종 상태를 추가로 저장하지 않았으므로 이후 실행도 단계 2에서 출발한다. 매 실행의 마지막 상태를 남기려면 반복문 뒤에 save를 호출한다. 그러면 다음 실행은 완료 단계 4를 복원하고 추가 계산 없이 요약만 출력한다.