<p><b>mpetersen@lanl.gov</b> 2011-09-30 14:38:46 -0600 (Fri, 30 Sep 2011)</p><p>Revised split_explicit_timestepping branch so that barotropic<br>
subcycling has a great deal of freedom in the ordering and number of u<br>
and ssh solves.<br>
<br>
The code now has the ability to solve consecutively for <br>
u, ssh, u, ssh<br>
in each barotropic subcycle.<br>
<br>
For previous mpas-ocean method (u, ssh, u) use these flags:<br>
config_btr_gam1_uWt1 = 0.5<br>
config_btr_gam2_SSHWt1 = 1.0<br>
config_btr_gam3_uWt2 = 1.0<br>
config_btr_solve_SSH2 = .false.<br>
<br>
For Griffies 2004 method (ssh, u, ssh), use:<br>
config_btr_gam1_uWt1 = 0.0<br>
config_btr_gam2_SSHWt1 = 0.2<br>
config_btr_gam3_uWt2 = 1.0<br>
config_btr_solve_SSH2 = .true.<br>
<br>
I have a bit-for-bit match between previous code (timestepping_mrp<br>
branch) and the first set of flags now.<br>
<br>
See Stage 2 in section 3.2 (eqn 3.31) in the updated<br>
branches/ocean_projects/doc/time_splitting_design/time_split_design.pdf<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/doc/time_splitting_design/time_split_design.pdf
===================================================================
(Binary files differ)
Modified: branches/ocean_projects/doc/time_splitting_design/time_split_design.tex
===================================================================
--- branches/ocean_projects/doc/time_splitting_design/time_split_design.tex 2011-09-30 19:38:22 UTC (rev 1048)
+++ branches/ocean_projects/doc/time_splitting_design/time_split_design.tex 2011-09-30 20:38:46 UTC (rev 1049)
@@ -276,18 +276,22 @@
\right)\\
\mbox{boundary update on }{\tilde \ubtr}_{n+j/J}\\
\zeta_{n+(j-1)/J}^{edge} = Interp(\zeta_{n+(j-1)/J}) \\
-{\tilde \ubtr}^*_{n+j/J} = \frac{1}{2}({\ubtr}_{n+j/J} + {\tilde \ubtr}_{n+j/J}) \\
-{\bf F}_j = {\tilde \ubtr}^*_{n+j/J}
+\tilde{\bf F}_j = \left((1-\gamma_1){\ubtr}_{n+(j-1)/J} + \gamma_1{\tilde \ubtr}_{n+j/J}\right)
\left(\zeta_{n+(j-1)/J}^{edge} + \sum_{k=1}^{N^{edge}} \Delta z_k \right)\\
-\zeta_{n+j/J} = \zeta_{n+(j-1)/J} + \frac{\Delta t}{J} \left( -
- </font>
<font color="red">abla \cdot {\bf F}_j \right) \\
-\mbox{boundary update on } \zeta_{n+j/J} \\
+\tilde\zeta_{n+j/J} = \zeta_{n+(j-1)/J} + \frac{\Delta t}{J} \left( -
+ </font>
<font color="red">abla \cdot \tilde{\bf F}_j \right) \\
+\mbox{boundary update on } \tilde\zeta_{n+j/J} \\
\mbox{compute }{\tilde \ubtr}_{n+j/J}^{\perp}\mbox{ from }{\tilde \ubtr}_{n+j/J} \\
\ubtr_{n+j/J} = \ubtr_{n+(j-1)/J} + \frac{\Delta t}{J} \left(
- f{\tilde \ubtr}_{n+j/J}^{\perp}
-- g</font>
<font color="blue">abla \zeta_{n+j/J} + {\overline {\bf G}_j}
+- g</font>
<font color="red">abla \left((1-\gamma_2)\zeta_{n+(j-1)/J} + \gamma_2\tilde\zeta_{n+j/J}\right) + {\overline {\bf G}_j}
\right),\\
-\mbox{boundary update on }\ubtr_{n+j/J}
+\mbox{boundary update on }\ubtr_{n+j/J}\\
+{\bf F}_j = \left((1-\gamma_3){\ubtr}_{n+(j-1)/J} + \gamma_3 \ubtr_{n+j/J}\right)
+\left(\zeta_{n+(j-1)/J}^{edge} + \sum_{k=1}^{N^{edge}} \Delta z_k \right)\\
+\zeta_{n+j/J} = \zeta_{n+(j-1)/J} + \frac{\Delta t}{J} \left( -
+ </font>
<font color="gray">abla \cdot {\bf F}_j \right) \\
+\mbox{boundary update on } \zeta_{n+j/J} \\
\end{array}
\right. \;\;\; j=1\ldots J \hspace{1cm}
\end{eqnarray}
Modified: branches/ocean_projects/split_explicit_timestepping/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/split_explicit_timestepping/src/core_ocean/mpas_ocn_time_integration_split.F 2011-09-30 19:38:22 UTC (rev 1048)
+++ branches/ocean_projects/split_explicit_timestepping/src/core_ocean/mpas_ocn_time_integration_split.F 2011-09-30 20:38:46 UTC (rev 1049)
@@ -87,7 +87,7 @@
type (block_type), pointer :: block
real (kind=RKIND) :: uhSum, hSum, sshEdge, flux, &
- uPerp, uCorr, tracerTemp, coef
+ uPerp, uCorr, tracerTemp, coef, FBtr_coeff, sshCell1, sshCell2
real (kind=RKIND), dimension(:), pointer :: sshNew
integer :: num_tracers, ucorr_coef, err
@@ -350,6 +350,7 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Barotropic subcycle: initial solve for velecity
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if (config_btr_gam1_uWt1>1.0e-12) then ! only do this part if it is needed in next SSH solve
uPerpTime = oldBtrSubcycleTime
block => domain % blocklist
@@ -423,19 +424,30 @@
block => block % next
end do ! block
+ endif ! config_btr_gam1_uWt1>1.0e-12
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Barotropic subcycle: Compute thickness flux and new SSH
+ ! Barotropic subcycle: Compute thickness flux and new SSH: PREDICTOR
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
block => domain % blocklist
do while (associated(block))
block % tend % ssh % array(:) = 0.0
- ! config_btr_flux_coef sets the forward weighting of velocity in the SSH computation
- ! config_btr_flux_coef= 1 flux = uBtrNew*H
- ! config_btr_flux_coef=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
- ! config_btr_flux_coef= 0 flux = uBtrOld*H
+ if (config_btr_solve_SSH2) then
+ ! If config_btr_solve_SSH2=.true., then do NOT accumulate FBtr in this SSH predictor
+ ! section, because it will be accumulated in the SSH corrector section.
+ FBtr_coeff = 0.0
+ else
+ ! otherwise, DO accumulate FBtr in this SSH predictor section
+ FBtr_coeff = 1.0
+ endif
+
+ ! config_btr_gam1_uWt1 sets the forward weighting of velocity in the SSH computation
+ ! config_btr_gam1_uWt1= 1 flux = uBtrNew*H
+ ! config_btr_gam1_uWt1=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
+ ! config_btr_gam1_uWt1= 0 flux = uBtrOld*H
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -445,9 +457,9 @@
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
- flux = ((1.0-config_btr_flux_coef) &
+ flux = ((1.0-config_btr_gam1_uWt1) &
* block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- + config_btr_flux_coef &
+ + config_btr_gam1_uWt1 &
* block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
* block % mesh % dvEdge % array(iEdge) &
* (sshEdge + hSum)
@@ -457,7 +469,7 @@
block % state % time_levs(1) % state % FBtr % array(iEdge) &
= block % state % time_levs(1) % state % FBtr % array(iEdge) &
- + flux
+ + FBtr_coeff*flux
end do
! SSHnew = SSHold + dt/J*(-div(Flux))
@@ -473,7 +485,7 @@
block => block % next
end do ! block
- ! boundary update on SSNnew
+ ! boundary update on SSHnew
block => domain % blocklist
do while (associated(block))
@@ -487,21 +499,7 @@
block => block % next
end do ! block
- block => domain % blocklist
- do while (associated(block))
- do iCell=1,block % mesh % nCells
-
- ! Accumulate SSH in running sum over the subcycles.
- block % state % time_levs(2) % state % ssh % array(iCell) &
- = block % state % time_levs(2) % state % ssh % array(iCell) &
- + block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell)
-
- end do
-
- block => block % next
- end do ! block
-
! mrp 110801 begin
! This whole section, bounded by 'mrp 110801', may be deleted later if it is found
! that barotropic del2 is not useful.
@@ -612,18 +610,42 @@
else
! uBtrNew = uBtrOld + dt*(-f*uBtroldPerp - g*grad(SSH) + G)
+! mrp 110926 change to weighted avg of SSH. revised to:
+
+ sshCell1 = &
+ (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
+
+ sshCell2 = &
+ (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
= block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- + dt/config_n_btr_subcycles *( &
+ + config_btr_gam3_uWt2*dt/config_n_btr_subcycles *( &
uPerp &
- gravity &
- *( block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
- - block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
+ *( sshCell2 &
+ - sshCell1 )&
/block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge) &
+ block % state % time_levs(1) % state % u_diffusionBtr % array(iEdge))
! added del2 diffusion to btr solve
+! mrp 110926 change to weighted avg of SSH. orig:
+! block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+! = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+! + dt/config_n_btr_subcycles *( &
+! uPerp &
+! - gravity &
+! *( 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(1) % state % u_diffusionBtr % array(iEdge))
+ ! added del2 diffusion to btr solve
+! mrp 110926 change to weighted avg of SSH. end
+
endif
end do
@@ -662,13 +684,100 @@
end do !do BtrCorIter=1,config_n_btr_cor_iter
+! mrp 110922 for extra corrector eta at end
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: Compute thickness flux and new SSH: CORRECTOR
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if (config_btr_solve_SSH2) then
+ block => domain % blocklist
+ do while (associated(block))
+
+ block % tend % ssh % array(:) = 0.0
+
+ ! config_btr_gam3_uWt2 sets the forward weighting of velocity in the SSH computation
+ ! config_btr_gam3_uWt2= 1 flux = uBtrNew*H
+ ! config_btr_gam3_uWt2=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
+ ! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
+
+ do iEdge=1,block % mesh % nEdges
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ sshEdge = 0.5 &
+ *( block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
+ hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
+
+ flux = ((1.0-config_btr_gam3_uWt2) &
+ * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam3_uWt2 &
+ * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * block % mesh % dvEdge % array(iEdge) &
+ * (sshEdge + hSum)
+
+ block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux
+ block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux
+
+ block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ = block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ + flux
+
+ end do
+
+ ! SSHnew = SSHold + dt/J*(-div(Flux))
+ do iCell=1,block % mesh % nCells
+
+ block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ + dt/config_n_btr_subcycles &
+ * block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
+
+ end do
+
+ block => block % next
+ end do ! block
+
+ ! boundary update on SSHnew
+ block => domain % blocklist
+ do while (associated(block))
+
+ call dmpar_exch_halo_field1dReal(domain % dminfo, &
+ block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
+ block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
+ block => block % next
+ end do ! block
+
+! mrp 110923 print
+!print *, 'sshsubcycleOld 2',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 *, 'sshSubcycleNew 2',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))
+! mrp 110923 print end
+
+ endif ! config_btr_solve_SSH2
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: Accumulate running sums, advance timestep pointers
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ ! Accumulate SSH in running sum over the subcycles.
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(2) % state % ssh % array(iCell) &
+ = block % state % time_levs(2) % state % ssh % array(iCell) &
+ + block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell)
+ end do
+
! 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) &
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-09-30 19:38:22 UTC (rev 1048)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-09-30 20:38:46 UTC (rev 1049)
@@ -186,9 +186,9 @@
var persistent real h ( nVertLevels nCells Time ) 2 ir h state - -
var persistent real rho ( nVertLevels nCells Time ) 2 ir rho state - -
var persistent real temperature ( nVertLevels nCells Time ) 2 iro temperature state tracers dynamics
-var persistent real salinity ( nVertLevels nCells Time ) 2 ir salinity state tracers dynamics
-var persistent real tracer1 ( nVertLevels nCells Time ) 2 ir tracer1 state tracers testing
-var persistent real tracer2 ( nVertLevels nCells Time ) 2 ir tracer2 state tracers testing
+var persistent real salinity ( nVertLevels nCells Time ) 2 iro salinity state tracers dynamics
+var persistent real tracer1 ( nVertLevels nCells Time ) 2 iro tracer1 state tracers testing
+var persistent real tracer2 ( nVertLevels nCells Time ) 2 iro tracer2 state tracers testing
# Tendency variables: neither read nor written to any files
var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
</font>
</pre>