<p><b>mpetersen@lanl.gov</b> 2011-06-02 16:24:55 -0600 (Thu, 02 Jun 2011)</p><p>Splitting is complete, and is tested for the following cases. I'm not sure if the split version works for other domains yet, and have not yet run it out in time.<br>
<br>
The code is still full of debuggins print statements, and comments need to be cleaned up yet.<br>
<br>
Testing splitting: progress so far.<br>
<br>
p30a test using a two layers, flat bottom, double-gyre domain, 10s steps<br>
/usr/projects/climate/mpeterse/coy/mpas_100407/branches/ocean_projects/grid_gen_ocean/pop_to_mpas/110214_double_gyre_15x30x2_200km_strat/grid.nc<br>
<br>
avg KE after ten steps, with:<br>
<br>
1.1089986020660202E-09 unsplit<br>
1.1025874035438945E-09 split, 1 big iter, 1 bcl subcycles, with u correction, 1 btr Cor iter<br>
1.1001634565912827E-09 split, 1 big iter, 10 bcl subcycles, with u correction, 1 btr Cor iter<br>
1.0999296352309274E-09 split, 1 big iter, 100 bcl subcycles, with u correction, 1 btr Cor iter<br>
1.1001617140468155E-09 split, 1 big iter, 10 bcl subcycles, with u correction, 10 btr Cor iter<br>
1.1089845937972990E-09 split, 1 big iter, 10 bcl subcycles, without u correction, 1 btr Cor iter<br>
1.0982773144895987E-09 split, 2 big iter, 10 bcl subcycles, with u correction, 1 btr Cor iter<br>
1.0982768375869806E-09 split, 10 big iter, 10 bcl subcycles, with u correction, 1 btr Cor iter<br>
<br>
I did not look at plots in time or cross-sections. However, this process helped me find several bugs,<br>
and I believe that the model is at least working reasonably in all these different configurations.<br>
<br>
Also, min/max of ssh computed by subcycles or one large step with accumulated flux match to all digits<br>
on these runs.<br>
<br>
p30b moved to 25 layer, flat bottom, double-gyre domain, 10s steps<br>
avg KE after ten steps, with:<br>
3.7077584769060971E-10 unsplit<br>
3.7121089556771133E-10 split, 1 big iter, 1 bcl subcycles, with u correction, 1 btr Cor iter<br>
3.7080783720286676E-10 split, 2 big iter, 10 bcl subcycles, with u correction, 2 btr Cor iter<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/Makefile
===================================================================
--- branches/ocean_projects/timesplitting_mrp/Makefile 2011-06-02 16:52:36 UTC (rev 871)
+++ branches/ocean_projects/timesplitting_mrp/Makefile 2011-06-02 22:24:55 UTC (rev 872)
@@ -86,7 +86,7 @@
"CC = gcc" \
"SFC = ifort" \
"SCC = gcc" \
- "FFLAGS = -real-size 64 -O0 -convert big_endian " \
+ "FFLAGS = -real-size 64 -O3 -convert big_endian" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
Modified: branches/ocean_projects/timesplitting_mrp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-06-02 16:52:36 UTC (rev 871)
+++ branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-06-02 22:24:55 UTC (rev 872)
@@ -22,10 +22,12 @@
/
&timestep_higdon
- config_n_ts_iter = 10
+ config_n_ts_iter = 2
config_n_bcl_iter_beg = 1
- config_n_bcl_iter_mid = 1
- config_n_bcl_iter_end = 1
+ config_n_bcl_iter_mid = 2
+ config_n_bcl_iter_end = 2
+ config_n_btr_subcycles = 10
+ config_n_btr_cor_iter = 2
config_compute_tr_midstage = .true.
/
&hmix
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-06-02 16:52:36 UTC (rev 871)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-06-02 22:24:55 UTC (rev 872)
@@ -21,6 +21,7 @@
namelist integer timestep_higdon config_n_bcl_iter_mid 4
namelist integer timestep_higdon config_n_bcl_iter_end 4
namelist integer timestep_higdon config_n_btr_subcycles 10
+namelist integer timestep_higdon config_n_btr_cor_iter 1
namelist logical timestep_higdon config_compute_tr_midstage true
namelist logical sw_model config_h_ScaleWithMesh false
namelist real hmix config_h_mom_eddy_visc2 0.0
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-06-02 16:52:36 UTC (rev 871)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-06-02 22:24:55 UTC (rev 872)
@@ -83,8 +83,9 @@
call compute_mesh_scaling(mesh)
call rbfInterp_initialize(mesh)
- call init_reconstruct(mesh)
- call reconstruct(block % state % time_levs(1) % state, mesh)
+! mrp temp
+! call init_reconstruct(mesh)
+! call reconstruct(block % state % time_levs(1) % state, mesh)
! initialize velocities and tracers on land to be -1e34
! The reconstructed velocity on land will have values not exactly
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-06-02 16:52:36 UTC (rev 871)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-06-02 22:24:55 UTC (rev 872)
@@ -654,7 +654,7 @@
= block % state % time_levs(1) % state % uBtr % array(iEdge)
! uBtrNew = BtrOld This is the first for the summation
- block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ block % state % time_levs(2) % state % uBtr % array(iEdge) &
= block % state % time_levs(1) % state % uBtr % array(iEdge)
! FBtr = 0
@@ -677,6 +677,8 @@
!! mrp del allocate(tend_ssh(block % mesh % nCells))
!print *, '7.1'
+ block % tend % ssh % array(:) = 0.0
+
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -755,20 +757,28 @@
! uBtrNew(iEdge) = uBtrNew(iEdge) + uBtrSubcycleNew(iEdge)
! end do
+ do BtrCorIter=1,config_n_btr_cor_iter
+
+ if (BtrCorIter==1) then
+ uPerpTime = oldBtrSubcycleTime
+ else
+ uPerpTime = newBtrSubcycleTime
+ endif
+
block => domain % blocklist
do while (associated(block))
! Iterate on u two times.
- do BtrCorIter=1,2
- if (BtrCorIter==1) uPerpTime = oldBtrSubcycleTime
- if (BtrCorIter==2) uPerpTime = newBtrSubcycleTime
+
+!print *, '9.1'
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+!print *, '9.2'
! Compute -f*uPerp
uPerp = 0.0
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
@@ -777,6 +787,7 @@
* block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
* block % mesh % fEdge % array(eoe)
end do
+!print *, '9.3'
! timestep in uBtrSubcycle:
! \ubtr_{n+(j+1)/J} = \ubtr_{n+j/J} + \frac{\Delta t}{J} \left(
@@ -785,24 +796,20 @@
! \right),
! uBtrSubCycleNew(iEdge) = uBtrSubcycleOld(iEdge) + dt*(fEdge(iEdge)*uBtrPerp(iEdge)
! - grad_ssh + GBtrForcing(iEdge))
+
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
= block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ dt/config_n_btr_subcycles *( &
uPerp &
- + block % state % time_levs(2) % state % u % array (k,iEdge) &
- grav_rho0Inv &
*( block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
- block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
/block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
-
- block % state % time_levs(2) % state % uBtr % array(iEdge) &
- = block % state % time_levs(2) % state % uBtr % array(iEdge) &
- + block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
+print *, '9.4'
end do
+print *, '9.5'
- enddo !do BtrCorIter=1,2
-
! printing:
nCells = block % mesh % nCells
nEdges = block % mesh % nEdges
@@ -813,8 +820,10 @@
block => block % next
end do ! block
-!print *, '10'
+print *, '10'
+
+
! boundary update on \ubtr_{n+(j+1)/J}
block => domain % blocklist
do while (associated(block))
@@ -827,10 +836,29 @@
block => block % next
end do ! block
+ end do !do BtrCorIter=1,2
+
+
+ ! uBtrNew = uBtrNew + uBtrSubcycleNEW
+ ! This accumulates the sum.
+ ! If the Barotropic Coriolis iteration is limited to one, this could
+ ! be merged with the above code.
+ block => domain % blocklist
+ do while (associated(block))
+ do iEdge=1,block % mesh % nEdges
+
+ block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ + block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
+
+ end do ! iEdge
+ block => block % next
+ end do ! block
+
! advance time pointers
oldBtrSubcycleTime = mod(oldBtrSubcycleTime,2)+1
newBtrSubcycleTime = mod(newBtrSubcycleTime,2)+1
-!print *, '11 oldBtrSubcycleTime, newBtrSubcycleTime',oldBtrSubcycleTime, newBtrSubcycleTime
+print *, '11 oldBtrSubcycleTime, newBtrSubcycleTime',oldBtrSubcycleTime, newBtrSubcycleTime
end do ! j=1,config_n_btr_subcycles
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -883,6 +911,7 @@
! allocate(tend_ssh(block % mesh % nCells))
+ block % tend % ssh % array(:) = 0.0
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -901,6 +930,10 @@
end do
! Now can compare sshSubcycleNEW (big step using summed fluxes) with
! sshSubcycleOLD (individual steps to get there)
+print *, 'ssh, by substeps',minval(block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve)), &
+ maxval(block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve))
+print *, 'ssh, by 1 step ',minval(block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve)), &
+ maxval(block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve))
! Correction velocity, u^{corr}:
!u^{corr} = \left( {\overline {\bf F}}
@@ -927,8 +960,12 @@
uhSum = uhSum + block % mesh % hZLevel % array(k) * uTemp(k)
hSum = hSum + block % mesh % hZLevel % array(k)
enddo
- uCorr = (block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum
+! mrp temp 110602 testing split: remove correction term.
+ uCorr = (block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum
+! uCorr = 0.0
+! mrp temp 110602 testing split end
+
! put u^{tr}, the velocity for tracer transport, in uNew
block % state % time_levs(2) % state % u % array(:,iEdge) = uTemp(:) + uCorr
@@ -971,7 +1008,7 @@
-!print *, '6'
+print *, '16'
!
! Stage 3: Tracer, density, pressure, vertical velocity prediction
@@ -990,15 +1027,15 @@
! printing:
nCells = block % mesh % nCells
-print *, 'h_tend ',minval(block % tend % h % array(1,1:nCells)),&
+print *, '1h_tend ',minval(block % tend % h % array(1,1:nCells)),&
maxval(block % tend % h % array(1,1:nCells))
-print *, 'tr_tend',minval(block % tend % tracers % array(3,1,1:nCells)),&
+print *, '1tr_tend',minval(block % tend % tracers % array(3,1,1:nCells)),&
maxval(block % tend % tracers % array(3,1,1:nCells))
block => block % next
end do
-!print *, '7'
+print *, '17'
! --- update halos for prognostic variables
@@ -1012,7 +1049,7 @@
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
-!print *, '8'
+print *, '18'
block => domain % blocklist
@@ -1296,7 +1333,8 @@
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
- call reconstruct(block % state % time_levs(2) % state, block % mesh)
+! mrp temp
+! call reconstruct(block % state % time_levs(2) % state, block % mesh)
block => block % next
end do
</font>
</pre>