From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 93003 invoked by alias); 19 Sep 2018 17:04:08 -0000 Mailing-List: contact fortran-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Post: List-Help: , Sender: fortran-owner@gcc.gnu.org Received: (qmail 91981 invoked by uid 89); 19 Sep 2018 17:04:07 -0000 Authentication-Results: sourceware.org; auth=none X-Spam-SWARE-Status: No, score=-4.1 required=5.0 tests=AWL,BAYES_05,GIT_PATCH_2,KAM_LAZY_DOMAIN_SECURITY,KAM_SHORT autolearn=ham version=3.3.2 spammy=West, East, departure, 3600 X-HELO: moene.org Received: from moene.org (HELO moene.org) (80.101.130.238) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with ESMTP; Wed, 19 Sep 2018 17:04:04 +0000 Received: from localhost ([127.0.0.1] helo=moene.org) by moene.org with esmtp (Exim 4.91) (envelope-from ) id 1g2ftQ-0001ph-HO; Wed, 19 Sep 2018 19:04:00 +0200 To: Thomas Koenig , Nicolas Koenig Cc: gfortran From: Toon Moene Subject: For your perusal when developing the new coarray code: A "mock" weather forecasting program using coarrays. Message-ID: Date: Wed, 19 Sep 2018 17:04:00 -0000 User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:52.0) Gecko/20100101 Thunderbird/52.9.1 MIME-Version: 1.0 Content-Type: multipart/mixed; boundary="------------E0E026FC54788D7C5DE46D47" X-IsSubscribed: yes X-SW-Source: 2018-09/txt/msg00133.txt.bz2 This is a multi-part message in MIME format. --------------E0E026FC54788D7C5DE46D47 Content-Type: text/plain; charset=utf-8; format=flowed Content-Transfer-Encoding: 7bit Content-length: 606 Thomas, Nicolas, The attached program doesn't do any physics, but it goes through the motions of computing and communications like the real thing. I tested it with opencoarrays, and it does a convincing run with the following command line (i.e., an empty configuration namelist): echo ' &config /' | cafrun -np 9 ./a.out Hope it is useful for you. -- Toon Moene - e-mail: toon@moene.org - phone: +31 346 214290 Saturnushof 14, 3738 XG Maartensdijk, The Netherlands At home: http://moene.org/~toon/; weather: http://moene.org/~hirlam/ Progress of GNU Fortran: http://gcc.gnu.org/wiki/GFortran#news --------------E0E026FC54788D7C5DE46D47 Content-Type: text/x-fortran; name="random-weather.f90" Content-Transfer-Encoding: 7bit Content-Disposition: attachment; filename="random-weather.f90" Content-length: 18283 ! ! Toon Moene 20180913-15,17-19 ! ! A mock "weather" forecasting program that goes through the motions ! of "representing" the weather on a limited area spread over several ! images to show the use of coarrays. ! module grid integer :: nxglobal[*], nyglobal[*], nzglobal[*], bdsz[*], flglobal[*], tsglobal[*] integer :: nxlocal, nylocal, boundary_size integer :: hgrid_step, vgrid_step, time_step end module grid ! module types ! ! Declaration of the type 'model_state' and its type bound procedures. ! type model_state real, allocatable :: ps(:,:)[:,:] ! Surface pressure real, allocatable :: t(:,:,:)[:,:] ! Temperature real, allocatable :: u(:,:,:)[:,:] ! U-component of wind (along parallels) real, allocatable :: v(:,:,:)[:,:] ! V-component of wind (along meridians) real, allocatable :: w(:,:,:)[:,:] ! W-component of wind (vertical) real, allocatable :: q(:,:,:)[:,:] ! Specific humidity contains procedure :: alloc ! Allocate a model state procedure :: dealloc ! Deallocate a model state procedure :: init ! Initialize a model state procedure :: random ! Construct a model state of random values -1..+1 procedure :: add ! Add two model states procedure :: ms_mult ! Multiply two model states procedure :: int_mult ! Multiply a model state with an integer procedure :: real_mult ! Multiply a model state with a real procedure :: assign ! Assign one model state to another procedure :: sla ! Semi-Lagrangian Advection Operator procedure :: bndrel ! Boundary relaxation end type model_state ! interface operator(+) module procedure add ! Standard end interface ! Overloading interface operator(*) ! Of module procedure ms_mult, & ! int_mult, real_mult ! Operators end interface ! And interface assignment(=) ! Assignment module procedure assign end interface interface operator(.sla.) ! Operator for module procedure sla ! Semi-Lagrangian end interface ! Advection ! contains ! ! Implementation of the module procedures on type 'model_state'. ! subroutine alloc(ms, nx, ny, nz, npx) class(model_state) :: ms integer, intent(in) :: nx, ny, nz, npx allocate(ms % ps(nx, ny)[npx, *]) allocate(ms % t(nx, ny, nz)[npx, *]) allocate(ms % u(nx, ny, nz)[npx, *]) allocate(ms % v(nx, ny, nz)[npx, *]) allocate(ms % w(nx, ny, nz)[npx, *]) allocate(ms % q(nx, ny, nz)[npx, *]) end subroutine alloc subroutine dealloc(ms) class(model_state) :: ms deallocate(ms % ps, ms % t, ms % u, ms % v, ms % w, ms % q) end subroutine dealloc subroutine init(ms, ps, t, u, v, w, q) class(model_state) :: ms real, intent(in) :: ps, t, u, v, w, q ms % ps = ps ms % t = t ms % u = u ms % v = v ms % w = w ms % q = q end subroutine init subroutine random(ms) ! Random numbers between -1 and +1 class(model_state) :: ms call random_number(ms % ps); ms % ps = 1.0 - 2.0 * ms % ps call random_number(ms % t ); ms % t = 1.0 - 2.0 * ms % t call random_number(ms % u ); ms % u = 1.0 - 2.0 * ms % u call random_number(ms % v ); ms % v = 1.0 - 2.0 * ms % v call random_number(ms % w ); ms % w = 1.0 - 2.0 * ms % w call random_number(ms % q ); ms % q = 1.0 - 2.0 * ms % q end subroutine random function add(ms1, ms2) class(model_state), allocatable :: add class(model_state), intent(in) :: ms1, ms2 integer :: nx, ny, nz, npx nx = size(ms1 % t, 1) ny = size(ms1 % t, 2) nz = size(ms1 % t, 3) npx = ucobound(ms1 % t, 1) - lcobound(ms1 % t, 1) + 1 ! Would be nice if it could be done this way: npx = cosize(ms1 % t, 1) allocate(add) call add % alloc(nx, ny, nz, npx) add % ps = ms1 % ps + ms2 % ps add % t = ms1 % t + ms2 % t add % u = ms1 % u + ms2 % u add % v = ms1 % v + ms2 % v add % w = ms1 % w + ms2 % w add % q = ms1 % q + ms2 % q end function add function ms_mult(ms1, ms2) class(model_state), allocatable :: ms_mult class(model_state), intent(in) :: ms1, ms2 integer :: nx, ny, nz, npx nx = size(ms1 % t, 1) ny = size(ms1 % t, 2) nz = size(ms1 % t, 3) npx = ucobound(ms1 % t, 1) - lcobound(ms1 % t, 1) + 1 allocate(ms_mult) call ms_mult % alloc(nx, ny, nz, npx) ms_mult % ps = ms1 % ps * ms2 % ps ms_mult % t = ms1 % t * ms2 % t ms_mult % u = ms1 % u * ms2 % u ms_mult % v = ms1 % v * ms2 % v ms_mult % w = ms1 % w * ms2 % w ms_mult % q = ms1 % q * ms2 % q end function ms_mult function int_mult(ms, ifactor) class(model_state), allocatable :: int_mult class(model_state), intent(in) :: ms integer, intent(in) :: ifactor integer :: nx, ny, nz, npx nx = size(ms % t, 1) ny = size(ms % t, 2) nz = size(ms % t, 3) npx = ucobound(ms % t, 1) - lcobound(ms % t, 1) + 1 allocate(int_mult) call int_mult % alloc(nx, ny, nz, npx) int_mult % ps = ms % ps * ifactor int_mult % t = ms % t * ifactor int_mult % u = ms % u * ifactor int_mult % v = ms % v * ifactor int_mult % w = ms % w * ifactor int_mult % q = ms % q * ifactor end function int_mult function real_mult(ms, rfactor) class(model_state), allocatable :: real_mult class(model_state), intent(in) :: ms real, intent(in) :: rfactor integer :: nx, ny, nz, npx nx = size(ms % t, 1) ny = size(ms % t, 2) nz = size(ms % t, 3) npx = ucobound(ms % t, 1) - lcobound(ms % t, 1) + 1 allocate(real_mult) call real_mult % alloc(nx, ny, nz, npx) real_mult % ps = ms % ps * rfactor real_mult % t = ms % t * rfactor real_mult % u = ms % u * rfactor real_mult % v = ms % v * rfactor real_mult % w = ms % w * rfactor real_mult % q = ms % q * rfactor end function real_mult subroutine assign(ms1, ms2) class(model_state), intent(out) :: ms1 ! Note: Components deallocated on entry. class(model_state), intent(in ) :: ms2 integer :: nx, ny, nz, npx nx = size(ms2 % t, 1) ny = size(ms2 % t, 2) nz = size(ms2 % t, 3) npx = ucobound(ms2 % t, 1) - lcobound(ms2 % t, 1) + 1 call ms1 % alloc(nx, ny, nz, npx) ms1 % ps = ms2 % ps ms1 % t = ms2 % t ms1 % u = ms2 % u ms1 % v = ms2 % v ms1 % w = ms2 % w ms1 % q = ms2 % q end subroutine assign function sla(ms) use grid class(model_state), allocatable :: sla class(model_state), intent(in) :: ms integer :: nx, jx, ny, jy, nz, jz, npx, ioff(3), idep(3), ico(2) nx = size(ms % t, 1) ny = size(ms % t, 2) nz = size(ms % t, 3) npx = ucobound(ms % t, 1) - lcobound(ms % t, 1) + 1 ico = this_image(ms % t) allocate(sla) call sla % alloc(nx, ny, nz, npx) ! ! Semi-Lagrangian advection looks upstream (i.e., up*wind* - at the departure point) ! at the atmospheric quantities to determine their values in the arrival point. ! The real thing is more complicated than just transporting the values, but it ! gives an idea of the computations and *communications* involved. ! do jz = 1, nz do jy = 1, ny do jx = 1, nx ioff(1) = nint(ms % u(jx, jy, jz) / hgrid_step * time_step) idep(1) = max(1, min(nxglobal, jx + nxlocal * (ico(1)-1) - ioff(1))) - 1 ! Zero based ! ioff(2) = nint(ms % v(jx, jy, jz) / hgrid_step * time_step) idep(2) = max(1, min(nyglobal, jy + nylocal * (ico(2)-1) - ioff(2))) - 1 ! Zero based ! ioff(3) = nint(ms % w(jx, jy, jz) / vgrid_step * time_step) idep(3) = max(1, min(nz, jz - ioff(3))) ! One based ! if (jz == 1) & sla % ps(jx, jy) = & ms % ps(mod(idep(1),nxlocal)+1, mod(idep(2),nylocal)+1)[idep(1)/nxlocal+1, idep(2)/nylocal+1] sla % t(jx, jy, jz) = & ms % t(mod(idep(1),nxlocal)+1, mod(idep(2),nylocal)+1, idep(3))[idep(1)/nxlocal+1, idep(2)/nylocal+1] sla % u(jx, jy, jz) = & ms % u(mod(idep(1),nxlocal)+1, mod(idep(2),nylocal)+1, idep(3))[idep(1)/nxlocal+1, idep(2)/nylocal+1] sla % v(jx, jy, jz) = & ms % v(mod(idep(1),nxlocal)+1, mod(idep(2),nylocal)+1, idep(3))[idep(1)/nxlocal+1, idep(2)/nylocal+1] sla % w(jx, jy, jz) = & ms % w(mod(idep(1),nxlocal)+1, mod(idep(2),nylocal)+1, idep(3))[idep(1)/nxlocal+1, idep(2)/nylocal+1] sla % q(jx, jy, jz) = & ms % q(mod(idep(1),nxlocal)+1, mod(idep(2),nylocal)+1, idep(3))[idep(1)/nxlocal+1, idep(2)/nylocal+1] enddo enddo enddo end function sla subroutine bndrel(ms, ps, t, u, v, w, q) use grid class(model_state) :: ms real, intent(in) :: ps, t, u, v, w, q real :: alpha integer :: nx, jx, ny, jy, nz, jz, ibegin, iend, ico(2), uco(2) nx = size(ms % t, 1) ny = size(ms % t, 2) nz = size(ms % t, 3) ico = this_image(ms % ps) uco = ucobound(ms % ps) ! ! Boundary relaxation is done to keep the "weather" inside the limited domain ! current with the developments outside it. In this "mock" weather forecasting ! program, we relax to fixed mean values. ! ! To prevent the relaxation to be doubly performed at the corners, ! the following convention is followed: ! ! North ! |------- ! West | | East ! -------| ! South ! ! The origin of the grid is in the South-West corner. ! ! 1. Relax (part of) the South boundary, if present in this slab if (ico(2) == 1) then if (ico(1) == uco(1)) then iend = nx - boundary_size else iend = nx endif do jz = 1, nz do jy = 1, boundary_size alpha = real(boundary_size - jy + 1) / boundary_size do jx = 1, iend if (jz == 1) & ms % ps(jx, jy) = alpha * ps + (1.0 - alpha) * ms % ps(jx, jy) ms % t(jx, jy, jz) = alpha * t + (1.0 - alpha) * ms % t(jx, jy, jz) ms % u(jx, jy, jz) = alpha * u + (1.0 - alpha) * ms % u(jx, jy, jz) ms % v(jx, jy, jz) = alpha * v + (1.0 - alpha) * ms % v(jx, jy, jz) ms % w(jx, jy, jz) = alpha * w + (1.0 - alpha) * ms % w(jx, jy, jz) ms % q(jx, jy, jz) = alpha * q + (1.0 - alpha) * ms % q(jx, jy, jz) enddo enddo enddo endif ! 2. Relax (part of) the North boundary, if present in this slab if (ico(2) == uco(2)) then if (ico(1) == 1) then ibegin = boundary_size + 1 else ibegin = 1 endif do jz = 1, nz do jy = ny, ny - boundary_size alpha = real(jy - ny + boundary_size) / boundary_size do jx = ibegin, nx if (jz == 1) & ms % ps(jx, jy) = alpha * ps + (1.0 - alpha) * ms % ps(jx, jy) ms % t(jx, jy, jz) = alpha * t + (1.0 - alpha) * ms % t(jx, jy, jz) ms % u(jx, jy, jz) = alpha * u + (1.0 - alpha) * ms % u(jx, jy, jz) ms % v(jx, jy, jz) = alpha * v + (1.0 - alpha) * ms % v(jx, jy, jz) ms % w(jx, jy, jz) = alpha * w + (1.0 - alpha) * ms % w(jx, jy, jz) ms % q(jx, jy, jz) = alpha * q + (1.0 - alpha) * ms % q(jx, jy, jz) enddo enddo enddo endif ! 3. Relax (part of) the West boundary, if present in this slab if (ico(1) == 1) then if (ico(2) == 1) then ibegin = boundary_size + 1 else ibegin = 1 endif do jz = 1, nz do jy = ibegin, ny do jx = 1, boundary_size alpha = real(boundary_size - jx + 1) / boundary_size if (jz == 1) & ms % ps(jx, jy) = alpha * ps + (1.0 - alpha) * ms % ps(jx, jy) ms % t(jx, jy, jz) = alpha * t + (1.0 - alpha) * ms % t(jx, jy, jz) ms % u(jx, jy, jz) = alpha * u + (1.0 - alpha) * ms % u(jx, jy, jz) ms % v(jx, jy, jz) = alpha * v + (1.0 - alpha) * ms % v(jx, jy, jz) ms % w(jx, jy, jz) = alpha * w + (1.0 - alpha) * ms % w(jx, jy, jz) ms % q(jx, jy, jz) = alpha * q + (1.0 - alpha) * ms % q(jx, jy, jz) enddo enddo enddo endif ! 4. Relax (part of) the East boundary, if present in this slab if (ico(1) == uco(1)) then if (ico(2) == uco(2)) then iend = ny - boundary_size else iend = ny endif do jz = 1, nz do jy = 1, iend do jx = nx, nx - boundary_size alpha = real(jx - nx + boundary_size) / boundary_size if (jz == 1) & ms % ps(jx, jy) = alpha * ps + (1.0 - alpha) * ms % ps(jx, jy) ms % t(jx, jy, jz) = alpha * t + (1.0 - alpha) * ms % t(jx, jy, jz) ms % u(jx, jy, jz) = alpha * u + (1.0 - alpha) * ms % u(jx, jy, jz) ms % v(jx, jy, jz) = alpha * v + (1.0 - alpha) * ms % v(jx, jy, jz) ms % w(jx, jy, jz) = alpha * w + (1.0 - alpha) * ms % w(jx, jy, jz) ms % q(jx, jy, jz) = alpha * q + (1.0 - alpha) * ms % q(jx, jy, jz) enddo enddo enddo endif end subroutine bndrel ! end module types ! program random_weather ! ! 0. Declarations. ! ! 0.1. Modules ! use grid use types ! implicit none ! ! 0.2. Parameters ! Pascal Kelvin m/s m/s m/s kg/kg real, parameter :: PS = 100000.0, T = 300.0, U = 0.0, V = 0.0, W = 0.0, Q = 0.002 ! Mean value real, parameter :: & dPSdt = 0.1, dTdt = 0.002, dUdt = 0.1, dVdt = 0.1, dWdt = 0.001, dQdt = 0.000001 ! Maximum tendency ! ! 0.3. Configuration variables ! seconds integer :: nx, ny, nz, forecast_length namelist /config/ nx, ny, nz, boundary_size, hgrid_step, vgrid_step, forecast_length, time_step ! ! 0.4. Global variables ! type(model_state) :: ms_full, ms_tend, ms_random ! ! 0.5. Local variables ! integer :: i, time, numimages, thisimage, npx, npy ! ! 1. Determine the number of images and compute the domain decomposition. ! ! 1.0. Set the domain defaults ! nx = 90; ny = 90; nz = 30; hgrid_step = 10000; vgrid_step = 100; boundary_size = 8; forecast_length = 3600; time_step = 240 ! ! 1.1. Read the domain definition and distribute it. ! numimages = num_images() thisimage = this_image() if (thisimage == 1) then read(*,config) ! Why won't this work as nxglobal[:] = nx ? do i = 1, numimages nxglobal[i] = nx; nyglobal[i] = ny; nzglobal[i] = nz; bdsz[i] = boundary_size flglobal[i] = forecast_length; tsglobal[i] = time_step enddo endif sync all if (thisimage /= 1) then nx = nxglobal; ny = nyglobal; nz = nzglobal; boundary_size = bdsz forecast_length = flglobal; time_step = tsglobal endif ! ! 1.2. Compute the domain decomposition. Note: only horizontal decomposition. ! Try a decomposition that's close to square (obviously, this algorithm ! has to be improved). ! npx = int(sqrt(real(numimages))) npy = numimages / npx; if (npx * npy < numimages) npx = npx + 1 nxlocal = nx / npx; if (nxlocal * npx < nx) nxlocal = nxlocal + 1 nylocal = ny / npy; if (nylocal * npy < ny) nylocal = nylocal + 1 print '(5(a,i4),a)','Decomposition information on image ',thisimage,' is ', & npx,' * ',npy,' slabs with ',nxlocal,' * ',nylocal, & ' grid cells on this image.' ! ! 1.3. Allocate the model states. ! call ms_full % alloc(nxlocal, nylocal, nz, npx) call ms_tend % alloc(nxlocal, nylocal, nz, npx) call ms_random % alloc(nxlocal, nylocal, nz, npx) ! ! 2. "Compute" the initial states of the model "atmosphere". ! ! 2.1. The initial values of the atmospheric quantities. ! call ms_full % init(PS, T, U, V, W, Q) ! ! 2.2. The maximum values of the tendencies of the atmospheric quantities. ! call ms_tend % init(dPSdt, dTdt, dUdt, dVdt, dWdt, dQdt) ! ! 3. "Integrate" the model in time. ! do time = 0, forecast_length, time_step ! ! 3.1. Compute the (random) factor on the tendencies of the atmospheric quantities ! "due" to sub grid scale processes. ! call ms_random % random ! ! 3.2. Perform "semi-lagrangian" advection (i.e., downstream) and add the increments. ! ms_full = .sla. ms_full + ms_random * ms_tend * time_step ! ! 3.3. Perform boundary relaxation towards mean values (to keep this simple). ! call ms_full % bndrel(PS, T, U, V, W, Q) ! ! 3.4. Print some diagnostic model output. ! print*, 'Time',time,' Image',thisimage, & ' PS=', ms_full % ps(nxlocal/2, nylocal/2), & ' T= ', ms_full % t (nxlocal/2, nylocal/2, 1), & ' U= ', ms_full % u (nxlocal/2, nylocal/2, 1), & ' V= ', ms_full % v (nxlocal/2, nylocal/2, 1), & ' W= ', ms_full % w (nxlocal/2, nylocal/2, 1), & ' Q= ', ms_full % q (nxlocal/2, nylocal/2, 1) ! end do ! ! 4. Complete the program. ! ! 4.1. Deallocate the model states. ! call ms_full % dealloc call ms_tend % dealloc call ms_random % dealloc ! end program random_weather --------------E0E026FC54788D7C5DE46D47--