<p><b>duda</b> 2010-07-13 17:38:14 -0600 (Tue, 13 Jul 2010)</p><p>BRANCH COMMIT<br>
<br>
Add initial set of changes to parallelize non-hydrostatic core:<br>
<br>
1) Add halo exchanges.<br>
2) In some cases, change loop bounds from nCellsSolve to nCells to avoid<br>
additional halo exchanges.<br>
3) Increase size of arrays allocated in time integration to include a dummy cell/edge/vertex <br>
at index nCells+1, nEdges+1, or nVertices+1.<br>
<br>
<br>
M src/core_nhyd_atmos/module_time_integration.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-07-13 21:28:52 UTC (rev 374)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-07-13 23:38:14 UTC (rev 375)
@@ -73,8 +73,8 @@
logical, parameter :: debug = .false.
! logical, parameter :: debug = .true.
logical, parameter :: debug_mass_conservation = .true.
- logical, parameter :: do_microphysics = .false.
- logical, parameter :: scalar_advection = .false.
+ logical, parameter :: do_microphysics = .true.
+ logical, parameter :: scalar_advection = .true.
real (kind=RKIND) :: domain_mass, scalar_mass, scalar_min, scalar_max
real (kind=RKIND) :: global_domain_mass, global_scalar_mass, global_scalar_min, global_scalar_max
@@ -100,6 +100,43 @@
block => domain % blocklist
do while (associated(block))
+! theta
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! scalars
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(1) % state % scalars % array(:,:,:), &
+ num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! pressure
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % pressure % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! vorticity
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+! divergence
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! pv_edge
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % pv_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+! rtheta_p
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_p % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! ru
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % ru % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ block => block % next
+ end do
+
+ block => domain % blocklist
+ do while (associated(block))
! We are setting values in the halo here, so no communications are needed.
! Alternatively, we could just set owned cells and edge values and communicate after this block loop.
call rk_integration_setup( block % time_levs(2) % state, block % time_levs(1) % state, block % mesh )
@@ -136,6 +173,22 @@
! we will need to communicate the momentum tendencies here - we want tendencies for all edges of owned cells
! because we are solving for all edges of owned cells
!***********************************
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % rho % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % w % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
block => domain % blocklist
do while (associated(block))
@@ -159,10 +212,30 @@
if(debug) write(0,*) ' acoustic step complete '
! will need communications here for rtheta_pp
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_pp % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rw_p % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % ru_p % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ block => block % next
+ end do
end do ! end of small stimestep loop
! will need communications here for rho_pp
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rho_pp % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
block => domain % blocklist
do while (associated(block))
@@ -197,6 +270,18 @@
block => block % next
end do
+ block => domain % blocklist
+ do while (associated(block))
+! For now, we do scalar halo updates later on...
+! call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(TEND) % scalars % array(:,:,:), &
+! num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % scalars % array(:,:,:), &
+! num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
else
write(0,*) ' no scalar advection '
@@ -214,7 +299,58 @@
! need communications here to fill out u, w, theta, p, and pp, scalars, etc
! so that they are available for next RK step or the first rk substep of the next timestep
+ block => domain % blocklist
+ do while (associated(block))
+! NB: if any of these cause differences, better to uncomment the version after qd_kessler?
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pressure % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % rho % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % rho_p % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % w % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_vertex % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_cell % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % ke % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % v % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % rho_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % scalars % array(:,:,:), &
+ num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
end do ! rk_step loop
! microphysics here...
@@ -227,6 +363,29 @@
end do
end if
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_p % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % exner % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pressure % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rt_diabatic_tend % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % scalars % array(:,:,:), &
+ num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
! if(debug) then
block => domain % blocklist
do while (associated(block))
@@ -506,8 +665,8 @@
real (kind=RKIND) :: smdiv, c2, rcv
real (kind=RKIND), dimension( grid % nVertLevels ) :: du
real (kind=RKIND), dimension( grid % nVertLevels + 1 ) :: dpzx
- real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells ) :: ts, rs
- real (kind=RKIND), dimension( grid % nVertLevels + 1 , grid % nCells ) :: ws
+ real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells+1 ) :: ts, rs
+ real (kind=RKIND), dimension( grid % nVertLevels + 1 , grid % nCells+1 ) :: ws
integer :: cell1, cell2, iEdge, iCell, k
real (kind=RKIND) :: pgrad, flux1, flux2, flux, resm, epssm
@@ -810,7 +969,7 @@
! recover owned-cell values in block
- if( iCell <= nCellsSolve ) then
+! if( iCell <= nCellsSolve ) then ! If using this test, more halo exchanges will be needed
if(debug) then
if( iCell == 479 ) then
@@ -855,7 +1014,7 @@
* (exner(k,iCell)-exner_base(k,iCell)))
end do
- end if
+! end if
end do
@@ -868,7 +1027,7 @@
cell1 = CellsOnEdge(1,iEdge)
cell2 = CellsOnEdge(2,iEdge)
- if( cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+! if( cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then ! If using this test, more halo exchanges will be needed
do k = 1, nVertLevels
ruAvg(k,iEdge) = ru(k,iEdge) + (ruAvg(k,iEdge) / float(ns))
@@ -888,7 +1047,7 @@
w(k,cell1) = w(k,cell1)+flux/AreaCell(cell1)
enddo
- end if
+! end if
enddo
@@ -1657,8 +1816,8 @@
! tendency for density
! divergence_ru may calculated in the diagnostic subroutine - it is temporary
- allocate(divergence_ru(nVertLevels, nCells))
- allocate(qtot(nVertLevels, nCells))
+ allocate(divergence_ru(nVertLevels, nCells+1))
+ allocate(qtot(nVertLevels, nCells+1))
divergence_ru(:,:) = 0.0
h_divergence(:,:) = 0.
@@ -1727,7 +1886,7 @@
! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
!
- allocate(rv(nVertLevels, nEdges))
+ allocate(rv(nVertLevels, nEdges+1))
rv(:,:) = 0.0
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
@@ -1820,10 +1979,10 @@
if ( h_mom_eddy_visc4 > 0.0 ) then
- allocate(delsq_divergence(nVertLevels, nCells))
- allocate(delsq_u(nVertLevels, nEdges))
- allocate(delsq_circulation(nVertLevels, nVertices))
- allocate(delsq_vorticity(nVertLevels, nVertices))
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
delsq_u(:,:) = 0.0
@@ -2101,7 +2260,7 @@
if ( h_mom_eddy_visc4 > 0.0 ) then
- allocate(delsq_theta(nVertLevels, nCells))
+ allocate(delsq_theta(nVertLevels, nCells+1))
delsq_theta(:,:) = 0.
@@ -2336,7 +2495,7 @@
if ( h_theta_eddy_visc4 > 0.0 ) then
- allocate(delsq_theta(nVertLevels, nCells))
+ allocate(delsq_theta(nVertLevels, nCells+1))
delsq_theta(:,:) = 0.
</font>
</pre>