* For your perusal when developing the new coarray code: A "mock" weather forecasting program using coarrays.
@ 2018-09-19 17:04 Toon Moene
0 siblings, 0 replies; only message in thread
From: Toon Moene @ 2018-09-19 17:04 UTC (permalink / raw)
To: Thomas Koenig, Nicolas Koenig; +Cc: gfortran
[-- Attachment #1: Type: text/plain, Size: 606 bytes --]
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
[-- Attachment #2: random-weather.f90 --]
[-- Type: text/x-fortran, Size: 18283 bytes --]
!
! 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
^ permalink raw reply [flat|nested] only message in thread
only message in thread, other threads:[~2018-09-19 17:04 UTC | newest]
Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2018-09-19 17:04 For your perusal when developing the new coarray code: A "mock" weather forecasting program using coarrays Toon Moene
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).