<p><b>dwj07@fsu.edu</b> 2012-09-18 11:04:21 -0600 (Tue, 18 Sep 2012)</p><p> <br>
-- BRANCH COMMIT --<br>
<br>
First cut at OpenMP in ocean core, using blocks.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_mpas_core.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_mpas_core.F 2012-09-18 16:37:51 UTC (rev 2162)
+++ branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_mpas_core.F 2012-09-18 17:04:21 UTC (rev 2163)
@@ -339,10 +339,15 @@
integer :: itimestep
real (kind=RKIND) :: dt
type (block_type), pointer :: block_ptr
+ type (block_type) :: block
type (MPAS_Time_Type) :: currTime
character(len=StrKIND) :: timeStamp
integer :: ierr
+
+ !$omp parallel default(shared)
+
+ !$omp single private(block_ptr)
! Eventually, dt should be domain specific
dt = config_dt
@@ -355,10 +360,15 @@
block_ptr => domain % blocklist
do while(associated(block_ptr))
+ block = block_ptr
+ !$omp task firstprivate(block)
call ocn_time_average_init(block_ptr % state % time_levs(1) % state)
+ !$omp end task
block_ptr => block_ptr % next
end do
+ !$omp taskwait
+
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
itimestep = 0
@@ -373,9 +383,13 @@
block_ptr => domain % blocklist
do while(associated(block_ptr))
+ block = block_ptr
+ !$omp task firstprivate(block)
call ocn_build_forcing_arrays(currTime, block_ptr % mesh, ierr)
+ !$omp end task
block_ptr => block_ptr % next
end do
+ !$omp taskwait
call mpas_timer_start("time integration", .false., timeIntTimer)
call mpas_timestep(domain, itimestep, dt, timeStamp)
@@ -384,9 +398,13 @@
! Move time level 2 fields back into time level 1 for next time step
block_ptr => domain % blocklist
do while(associated(block_ptr))
+ block = block_ptr
+ !$omp task firstprivate(block)
call mpas_shift_time_levels_state(block_ptr % state)
+ !$omp end task
block_ptr => block_ptr % next
end do
+ !$omp taskwait
if (mpas_is_alarm_ringing(clock, outputAlarmID, ierr=ierr)) then
call mpas_reset_clock_alarm(clock, outputAlarmID, ierr=ierr)
@@ -398,17 +416,25 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
+ block = block_ptr
+ !$omp task firstprivate(block)
call ocn_time_average_normalize(block_ptr % state % time_levs(1) % state)
+ !$omp end task
block_ptr => block_ptr % next
end do
+ !$omp taskwait
call ocn_write_output_frame(output_obj, output_frame, domain)
block_ptr => domain % blocklist
do while (associated(block_ptr))
+ block = block_ptr
+ !$omp task firstprivate(block)
call ocn_time_average_init(block_ptr % state % time_levs(1) % state)
+ !$omp end task
block_ptr => block_ptr % next
end do
+ !$omp taskwait
end if
if (mpas_is_alarm_ringing(clock, restartAlarmID, ierr=ierr)) then
@@ -422,6 +448,10 @@
end do
+ !$omp end single
+
+ !$omp end parallel
+
end subroutine mpas_core_run!}}}
subroutine ocn_write_output_frame(output_obj, output_frame, domain)!{{{
Modified: branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_time_integration.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_time_integration.F 2012-09-18 16:37:51 UTC (rev 2162)
+++ branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_time_integration.F 2012-09-18 17:04:21 UTC (rev 2163)
@@ -88,6 +88,8 @@
type (dm_info) :: dminfo
type (block_type), pointer :: block
+ real (kind=RKIND) :: nanCheck
+
if (rk4On) then
call ocn_time_integrator_rk4(domain, dt)
elseif (splitOn) then
@@ -98,7 +100,9 @@
do while (associated(block))
block % state % time_levs(2) % state % xtime % scalar = timeStamp
- if (isNaN(sum(block % state % time_levs(2) % state % u % array))) then
+ nanCheck = sum(block % state % time_levs(2) % state % u % array)
+
+ if(nanCheck /= nanCheck) then
write(0,*) 'Abort: NaN detected'
call mpas_dmpar_abort(dminfo)
endif
Modified: branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_time_integration_rk4.F 2012-09-18 16:37:51 UTC (rev 2162)
+++ branches/omp_blocks/openmp_test/src/core_ocean_blocks/mpas_ocn_time_integration_rk4.F 2012-09-18 17:04:21 UTC (rev 2163)
@@ -78,6 +78,7 @@
integer :: iCell, k, i, err
type (block_type), pointer :: block
+ type (block_type) :: block_d
integer :: rk_step
@@ -103,19 +104,23 @@
!
block => domain % blocklist
do while (associated(block))
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
- do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % maxLevelCell % array(iCell)
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % h % array(k,iCell)
+ block_d = block
+ !$omp task firstprivate(block_d)
+ block_d % state % time_levs(2) % state % u % array(:,:) = block_d % state % time_levs(1) % state % u % array(:,:)
+ block_d % state % time_levs(2) % state % h % array(:,:) = block_d % state % time_levs(1) % state % h % array(:,:)
+ do iCell=1,block_d % mesh % nCells ! couple tracers to h
+ do k=1,block_d % mesh % maxLevelCell % array(iCell)
+ block_d % state % time_levs(2) % state % tracers % array(:,k,iCell) = block_d % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ * block_d % state % time_levs(1) % state % h % array(k,iCell)
end do
end do
- call mpas_copy_state(block % provis, block % state % time_levs(1) % state)
+ call mpas_copy_state(block_d % provis, block_d % state % time_levs(1) % state)
+ !$omp end task
block => block % next
end do
+ !$omp taskwait
rk_weights(1) = dt/6.
rk_weights(2) = dt/3.
@@ -148,25 +153,29 @@
call mpas_timer_start("RK4-tendency computations")
block => domain % blocklist
do while (associated(block))
+ block_d = block
+ !$omp task firstprivate(block_d)
! mrp 111206 put ocn_wtop call at top for ALE
- call ocn_wtop(block % provis, block % provis, block % mesh)
+ call ocn_wtop(block_d % provis, block_d % provis, block_d % mesh)
if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, err)
+ call ocn_vmix_coefs(block_d % mesh, block_d % provis, block_d % diagnostics, err)
end if
- call ocn_tend_h(block % tend, block % provis, block % mesh)
- call ocn_tend_u(block % tend, block % provis, block % diagnostics, block % mesh)
+ call ocn_tend_h(block_d % tend, block_d % provis, block_d % mesh)
+ call ocn_tend_u(block_d % tend, block_d % provis, block_d % diagnostics, block_d % mesh)
! mrp 110718 filter btr mode out of u_tend
! still got h perturbations with just this alone. Try to set uBtr=0 after full u computation
if (config_rk_filter_btr_mode) then
- call ocn_filter_btr_mode_tend_u(block % tend, block % provis, block % mesh)
+ call ocn_filter_btr_mode_tend_u(block_d % tend, block_d % provis, block_d % mesh)
endif
- call ocn_tend_scalar(block % tend, block % provis, block % diagnostics, block % mesh, dt)
+ call ocn_tend_scalar(block_d % tend, block_d % provis, block_d % diagnostics, block_d % mesh, dt)
+ !$omp end task
block => block % next
end do
+ !$omp taskwait
call mpas_timer_stop("RK4-tendency computations")
! --- update halos for prognostic variables
@@ -183,42 +192,46 @@
if (rk_step < 4) then
block => domain % blocklist
do while (associated(block))
+ block_d = block
- block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
+ !$omp task firstprivate(block_d)
+ block_d % provis % u % array(:,:) = block_d % state % time_levs(1) % state % u % array(:,:) &
+ + rk_substep_weights(rk_step) * block_d % tend % u % array(:,:)
- block % provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % maxLevelCell % array(iCell)
- block % provis % tracers % array(:,k,iCell) = ( block % state % time_levs(1) % state % h % array(k,iCell) * &
- block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
- ) / block % provis % h % array(k,iCell)
+ block_d % provis % h % array(:,:) = block_d % state % time_levs(1) % state % h % array(:,:) &
+ + rk_substep_weights(rk_step) * block_d % tend % h % array(:,:)
+ do iCell=1,block_d % mesh % nCells
+ do k=1,block_d % mesh % maxLevelCell % array(iCell)
+ block_d % provis % tracers % array(:,k,iCell) = ( block_d % state % time_levs(1) % state % h % array(k,iCell) * &
+ block_d % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ + rk_substep_weights(rk_step) * block_d % tend % tracers % array(:,k,iCell) &
+ ) / block_d % provis % h % array(k,iCell)
end do
end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ block_d % provis % u % array(:,:) = block_d % state % time_levs(1) % state % u % array(:,:)
end if
if (config_prescribe_velocity) then
- block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ block_d % provis % u % array(:,:) = block_d % state % time_levs(1) % state % u % array(:,:)
end if
if (config_prescribe_thickness) then
- block % provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ block_d % provis % h % array(:,:) = block_d % state % time_levs(1) % state % h % array(:,:)
end if
- call ocn_diagnostic_solve(dt, block % provis, block % mesh)
+ call ocn_diagnostic_solve(dt, block_d % provis, block_d % mesh)
! Compute velocity transport, used in advection terms of h and tracer tendency
- block % provis % uTransport % array(:,:) &
- = block % provis % u % array(:,:) &
- + block % provis % uBolusGM % array(:,:)
+ block_d % provis % uTransport % array(:,:) &
+ = block_d % provis % u % array(:,:) &
+ + block_d % provis % uBolusGM % array(:,:)
+ !$omp end task
block => block % next
end do
+ !$omp taskwait
end if
call mpas_timer_stop("RK4-update diagnostic variables")
@@ -227,23 +240,27 @@
call mpas_timer_start("RK4-RK4 accumulate update")
block => domain % blocklist
do while (associated(block))
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
- + rk_weights(rk_step) * block % tend % u % array(:,:)
+ block_d = block
+ !$omp task firstprivate(block_d)
+ block_d % state % time_levs(2) % state % u % array(:,:) = block_d % state % time_levs(2) % state % u % array(:,:) &
+ + rk_weights(rk_step) * block_d % tend % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
- + rk_weights(rk_step) * block % tend % h % array(:,:)
+ block_d % state % time_levs(2) % state % h % array(:,:) = block_d % state % time_levs(2) % state % h % array(:,:) &
+ + rk_weights(rk_step) * block_d % tend % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % maxLevelCell % array(iCell)
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- ( block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell) )&
- / block % state % time_levs(2) % state % h % array(k, iCell)
+ do iCell=1,block_d % mesh % nCells
+ do k=1,block_d % mesh % maxLevelCell % array(iCell)
+ block_d % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ ( block_d % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ + rk_weights(rk_step) * block_d % tend % tracers % array(:,k,iCell) )&
+ / block_d % state % time_levs(2) % state % h % array(k, iCell)
end do
end do
+ !$omp end task
block => block % next
end do
+ !$omp taskwait
call mpas_timer_stop("RK4-RK4 accumulate update")
end do
@@ -260,9 +277,13 @@
call mpas_timer_start("RK4-implicit vert mix")
block => domain % blocklist
do while(associated(block))
+ block_d = block
+ !$omp task firstprivate(block_d)
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ !$omp end task
block => block % next
end do
+ !$omp taskwait
! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
! this leads to lack of volume conservation. It is required because halo updates in RK4 are only
@@ -278,6 +299,8 @@
block => domain % blocklist
do while (associated(block))
+ block_d = block
+ !$omp task firstprivate(block_d)
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
@@ -316,9 +339,11 @@
!TDR
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
+ !$omp end task
block => block % next
end do
+ !$omp taskwait
call mpas_timer_stop("RK4-cleaup phase")
call mpas_deallocate_provis_states(domain % blocklist)
</font>
</pre>