<p><b>duda</b> 2010-02-24 14:47:50 -0700 (Wed, 24 Feb 2010)</p><p>Add positive definite and monotonic limiters for scalar transport.<br>
Limiting is controlled by new namelist parameters config_positive_definite<br>
and config_monotonic.<br>
<br>
Note: the new code still needs to be modified to work with <br>
multiple blocks per process.<br>
<br>
<br>
M src/module_time_integration.F<br>
M Registry/Registry<br>
M namelist.input<br>
</p><hr noshade><pre><font color="gray">Modified: branches/hyd_model/Registry/Registry
===================================================================
--- branches/hyd_model/Registry/Registry 2010-02-24 21:06:01 UTC (rev 119)
+++ branches/hyd_model/Registry/Registry 2010-02-24 21:47:50 UTC (rev 120)
@@ -15,6 +15,8 @@
namelist integer sw_model config_number_of_sub_steps 4
namelist integer sw_model config_theta_adv_order 2
namelist integer sw_model config_scalar_adv_order 2
+namelist logical sw_model config_positive_definite false
+namelist logical sw_model config_monotonic true
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
Modified: branches/hyd_model/namelist.input
===================================================================
--- branches/hyd_model/namelist.input 2010-02-24 21:06:01 UTC (rev 119)
+++ branches/hyd_model/namelist.input 2010-02-24 21:47:50 UTC (rev 120)
@@ -11,8 +11,10 @@
config_h_theta_eddy_visc2 = 0.0
config_h_theta_eddy_visc4 = 0.0
config_v_theta_eddy_visc2 = 0.0
- config_theta_adv_order = 2
- config_scalar_adv_order = 2
+ config_theta_adv_order = 3
+ config_scalar_adv_order = 3
+ config_positive_definite = .false.
+ config_monotonic = .true.
/
&io
Modified: branches/hyd_model/src/module_time_integration.F
===================================================================
--- branches/hyd_model/src/module_time_integration.F 2010-02-24 21:06:01 UTC (rev 119)
+++ branches/hyd_model/src/module_time_integration.F 2010-02-24 21:47:50 UTC (rev 120)
@@ -262,9 +262,21 @@
block => domain % blocklist
do while (associated(block))
- call advance_scalars( block % intermediate_step(TEND), &
- block % time_levs(1) % state, block % time_levs(2) % state, &
- block % mesh, rk_timestep(rk_step) )
+ !
+ ! Note: The advance_scalars_mono routine can be used without limiting, and thus, encompasses
+ ! the functionality of the advance_scalars routine; however, it is noticeably slower,
+ ! so we keep the advance_scalars routine as well
+ !
+ if (rk_step < 3 .or. (.not. config_monotonic .and. .not. config_positive_definite)) then
+ call advance_scalars( block % intermediate_step(TEND), &
+ block % time_levs(1) % state, block % time_levs(2) % state, &
+ block % mesh, rk_timestep(rk_step) )
+ else
+ call advance_scalars_mono( block % intermediate_step(TEND), &
+ block % time_levs(1) % state, block % time_levs(2) % state, &
+ block % mesh, rk_timestep(rk_step), rk_step, 3, &
+ domain % dminfo, block % parinfo % cellsToSend, block % parinfo % cellsToRecv )
+ end if
block => block % next
end do
@@ -1313,7 +1325,12 @@
integer :: nVertLevels
real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw
+ real (kind=RKIND) :: coef_3rd_order
+ coef_3rd_order = 0.
+ if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
+ if (config_scalar_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
scalar_old => s_old % tracers % array
scalar_new => s_new % tracers % array
deriv_two => grid % deriv_two % array
@@ -1384,12 +1401,25 @@
if (uhAvg(k,iEdge) > 0) then
flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
else
flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
end if
+
+! old version of the above code, with coef_3rd_order assumed to be 1.0
+! if (uhAvg(k,iEdge) > 0) then
+! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
+! 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+! -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
+! else
+! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
+! 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+! -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
+! end if
scalar_tend(iTracer,k,cell1) = scalar_tend(iTracer,k,cell1) - flux/areaCell(cell1)
scalar_tend(iTracer,k,cell2) = scalar_tend(iTracer,k,cell2) + flux/areaCell(cell2)
@@ -1462,6 +1492,349 @@
end subroutine advance_scalars
+ subroutine advance_scalars_mono( tend, s_old, s_new, grid, dt, rk_step, rk_order, dminfo, cellsToSend, cellsToRecv)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed scalar tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_state), intent(in) :: tend
+ type (grid_state), intent(in) :: s_old
+ type (grid_state), intent(out) :: s_new
+ type (grid_meta), intent(in) :: grid
+ integer, intent(in) :: rk_step, rk_order
+ real (kind=RKIND), intent(in) :: dt
+ type (dm_info), intent(in) :: dminfo
+ type (exchange_list), pointer :: cellsToSend, cellsToRecv
+
+ integer :: i, iCell, iEdge, k, iTracer, cell_upwind, cell1, cell2
+ real (kind=RKIND) :: flux, tracer_edge, d2fdx2_cell1, d2fdx2_cell2
+ real (kind=RKIND) :: fdir, flux_upwind, h_flux_upwind, s_upwind
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalar_old, scalar_new, scalar_tend
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND), dimension(:,:), pointer :: uhAvg, h_old, h_new, wwAvg
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ real (kind=RKIND), dimension( grid % nTracers, grid % nEdges) :: h_flux
+ real (kind=RKIND), dimension( grid % nTracers, grid % nCells, 2 ) :: v_flux, v_flux_upwind, s_update
+ real (kind=RKIND), dimension( grid % nTracers, grid % nCells, 2 ) :: scale_out, scale_in
+ real (kind=RKIND), dimension( grid % nTracers ) :: s_max, s_min, s_max_update, s_min_update
+
+ integer :: nVertLevels, km0, km1, ktmp, kcp1, kcm1
+
+ real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw
+ real (kind=RKIND), parameter :: eps=1.e-20
+ real (kind=RKIND) :: coef_3rd_order
+
+ scalar_old => s_old % tracers % array
+ scalar_new => s_new % tracers % array
+ deriv_two => grid % deriv_two % array
+ uhAvg => grid % uhAvg % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ scalar_tend => tend % tracers % array
+ h_old => s_old % h % array
+ h_new => s_new % h % array
+ wwAvg => grid % wwAvg % array
+ areaCell => grid % areaCell % array
+
+ fnm => grid % fnm % array
+ fnp => grid % fnp % array
+ rdnw => grid % rdnw % array
+
+ nVertLevels = grid % nVertLevels
+
+ scalar_tend = 0. ! testing purposes - we have no sources or sinks
+
+ !
+ ! Runge Kutta integration, so we compute fluxes from scalar_new values, update starts from scalar_old
+ !
+
+ km1 = 1
+ km0 = 2
+ v_flux(:,:,km1) = 0.
+ v_flux_upwind(:,:,km1) = 0.
+ scale_out(:,:,:) = 1.
+ scale_in(:,:,:) = 1.
+
+ coef_3rd_order = 0.
+ if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
+ if (config_scalar_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
+ do k = 1, grid % nVertLevels
+ kcp1 = min(k+1,grid % nVertLevels)
+ kcm1 = max(k-1,1)
+
+! vertical flux
+
+ do iCell=1,grid % nCells
+
+ if (k < grid % nVertLevels) then
+ cell_upwind = k
+ if (wwAvg(k+1,iCell) >= 0) cell_upwind = k+1
+ do iTracer=1,grid % nTracers
+ v_flux(iTracer,iCell,km0) = dt * wwAvg(k+1,iCell) * &
+ (fnm(k+1) * scalar_new(iTracer,k+1,iCell) + fnp(k+1) * scalar_new(iTracer,k,iCell))
+ v_flux_upwind(iTracer,iCell,km0) = dt * wwAvg(k+1,iCell) * scalar_old(iTracer,cell_upwind,iCell)
+ v_flux(iTracer,iCell,km0) = v_flux(iTracer,iCell,km0) - v_flux_upwind(iTracer,iCell,km0)
+! v_flux(iTracer,iCell,km0) = 0. ! use only upwind - for testing
+ s_update(iTracer,iCell,km0) = scalar_old(iTracer,k,iCell) * h_old(k,iCell) &
+ - rdnw(k) * (v_flux_upwind(iTracer,iCell,km0) - v_flux_upwind(iTracer,iCell,km1))
+ end do
+ else
+ do iTracer=1,grid % nTracers
+ v_flux(iTracer,iCell,km0) = 0.
+ v_flux_upwind(iTracer,iCell,km0) = 0.
+ s_update(iTracer,iCell,km0) = scalar_old(iTracer,k,iCell) * h_old(k,iCell) &
+ - rdnw(k) * (v_flux_upwind(iTracer,iCell,km0) - v_flux_upwind(iTracer,iCell,km1))
+ end do
+ end if
+
+ end do
+
+! horizontal flux
+
+ if (config_scalar_adv_order == 2) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+ cell_upwind = cell2
+ if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
+ do iTracer=1,grid % nTracers
+ tracer_edge = 0.5 * (scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2))
+ h_flux(iTracer,iEdge) = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * tracer_edge
+ h_flux_upwind = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * scalar_old(iTracer,k,cell_upwind)
+ h_flux(iTracer,iEdge) = h_flux(iTracer,iEdge) - h_flux_upwind
+! h_flux(iTracer,iEdge) = 0. ! use only upwind - for testing
+ s_update(iTracer,cell1,km0) = s_update(iTracer,cell1,km0) - h_flux_upwind / grid % areaCell % array(cell1)
+ s_update(iTracer,cell2,km0) = s_update(iTracer,cell2,km0) + h_flux_upwind / grid % areaCell % array(cell2)
+ end do
+ end if
+ end do
+
+ else if (config_scalar_adv_order >= 3) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+ cell_upwind = cell2
+ if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
+ do iTracer=1,grid % nTracers
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iTracer,k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ if (uhAvg(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
+ 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ else
+ flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
+ 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ end if
+
+ h_flux(iTracer,iEdge) = dt * flux
+ h_flux_upwind = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * scalar_old(iTracer,k,cell_upwind)
+ h_flux(iTracer,iEdge) = h_flux(iTracer,iEdge) - h_flux_upwind
+! h_flux(iTracer,iEdge) = 0. ! use only upwind - for testing
+ s_update(iTracer,cell1,km0) = s_update(iTracer,cell1,km0) - h_flux_upwind / grid % areaCell % array(cell1)
+ s_update(iTracer,cell2,km0) = s_update(iTracer,cell2,km0) + h_flux_upwind / grid % areaCell % array(cell2)
+ end do
+ end if
+ end do
+
+ end if
+
+
+ if ( (rk_step == rk_order) .and. (config_monotonic .or. config_positive_definite) ) then
+
+!*************************************************************************************************************
+!--- limiter - we limit horizontal and vertical fluxes on level k
+!--- (these are h fluxes contributing to level k scalars, and v flux contributing to level k, k-1 tracers)
+
+ do iCell=1,grid % nCells
+
+ do iTracer=1,grid % nTracers
+
+ s_max(iTracer) = max(scalar_old(iTracer,k,iCell), scalar_old(iTracer,kcp1,iCell), scalar_old(iTracer,kcm1,iCell))
+ s_min(iTracer) = min(scalar_old(iTracer,k,iCell), scalar_old(iTracer,kcp1,iCell), scalar_old(iTracer,kcm1,iCell))
+ s_max_update(iTracer) = s_update(iTracer,iCell,km0)
+ s_min_update(iTracer) = s_update(iTracer,iCell,km0)
+
+ ! add in vertical flux to get max and min estimate
+ s_max_update(iTracer) = s_max_update(iTracer) &
+ - rdnw(k) * (max(0.,v_flux(iTracer,iCell,km0)) - min(0.,v_flux(iTracer,iCell,km1)))
+ s_min_update(iTracer) = s_min_update(iTracer) &
+ - rdnw(k) * (min(0.,v_flux(iTracer,iCell,km0)) - max(0.,v_flux(iTracer,iCell,km1)))
+
+ end do
+
+ do i = 1, grid % nEdgesOnCell % array(iCell) ! go around the edges of each cell
+ if (grid % cellsOnCell % array(i,iCell) > 0) then
+ do iTracer=1,grid % nTracers
+
+ s_max(iTracer) = max(scalar_old(iTracer,k,grid % cellsOnCell % array(i,iCell)), s_max(iTracer))
+ s_min(iTracer) = min(scalar_old(iTracer,k,grid % cellsOnCell % array(i,iCell)), s_min(iTracer))
+
+ iEdge = grid % EdgesOnCell % array (i,iCell)
+ if (iCell == cellsOnEdge(1,iEdge)) then
+ fdir = 1.0
+ else
+ fdir = -1.0
+ end if
+ flux = -fdir * h_flux(iTracer,iEdge)/grid % areaCell % array(iCell)
+ s_max_update(iTracer) = s_max_update(iTracer) + max(0.,flux)
+ s_min_update(iTracer) = s_min_update(iTracer) + min(0.,flux)
+
+ end do
+ end if
+
+ end do
+
+ if( config_positive_definite ) s_min(:) = 0.
+
+ do iTracer=1,grid % nTracers
+ scale_out (iTracer,iCell,km0) = 1.
+ scale_in (iTracer,iCell,km0) = 1.
+ s_max_update (iTracer) = s_max_update (iTracer) / h_new (k,iCell)
+ s_min_update (iTracer) = s_min_update (iTracer) / h_new (k,iCell)
+ s_upwind = s_update(iTracer,iCell,km0) / h_new(k,iCell)
+ if ( s_max_update(iTracer) > s_max(iTracer) .and. config_monotonic) &
+ scale_in (iTracer,iCell,km0) = max(0.,(s_max(iTracer)-s_upwind)/(s_max_update(iTracer)-s_upwind+eps))
+ if ( s_min_update(iTracer) < s_min(iTracer) ) &
+ scale_out (iTracer,iCell,km0) = max(0.,(s_upwind-s_min(iTracer))/(s_upwind-s_min_update(iTracer)+eps))
+ end do
+
+ end do ! end loop over cells to compute scale factor
+
+
+ call dmpar_exch_halo_field2dReal(dminfo, scale_out(:,:,1), &
+ grid % nTracers, grid % nCells, &
+ cellsToSend, cellsToRecv)
+ call dmpar_exch_halo_field2dReal(dminfo, scale_out(:,:,2), &
+ grid % nTracers, grid % nCells, &
+ cellsToSend, cellsToRecv)
+ call dmpar_exch_halo_field2dReal(dminfo, scale_in(:,:,1), &
+ grid % nTracers, grid % nCells, &
+ cellsToSend, cellsToRecv)
+ call dmpar_exch_halo_field2dReal(dminfo, scale_in(:,:,2), &
+ grid % nTracers, grid % nCells, &
+ cellsToSend, cellsToRecv)
+
+ ! rescale the horizontal fluxes
+
+ do iEdge = 1, grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+ do iTracer=1,grid % nTracers
+ flux = h_flux(iTracer,iEdge)
+ if (flux > 0) then
+ flux = flux * min(scale_out(iTracer,cell1,km0), scale_in(iTracer,cell2,km0))
+ else
+ flux = flux * min(scale_in(iTracer,cell1,km0), scale_out(iTracer,cell2,km0))
+ end if
+ h_flux(iTracer,iEdge) = flux
+ end do
+ end if
+ end do
+
+ ! rescale the vertical flux
+
+ do iCell=1,grid % nCells
+ do iTracer=1,grid % nTracers
+ flux = v_flux(iTracer,iCell,km1)
+ if (flux > 0) then
+ flux = flux * min(scale_out(iTracer,iCell,km0), scale_in(iTracer,iCell,km1))
+ else
+ flux = flux * min(scale_in(iTracer,iCell,km0), scale_out(iTracer,iCell,km1))
+ end if
+ v_flux(iTracer,iCell,km1) = flux
+ end do
+ end do
+
+! end of limiter
+!*******************************************************************************************************************
+
+ end if
+
+!--- update
+
+ do iCell=1,grid % nCells
+ ! add in upper vertical flux that was just renormalized
+ do iTracer=1,grid % nTracers
+ s_update(iTracer,iCell,km0) = s_update(iTracer,iCell,km0) + rdnw(k) * v_flux(iTracer,iCell,km1)
+ if (k > 1) s_update(iTracer,iCell,km1) = s_update(iTracer,iCell,km1) - rdnw(k-1)*v_flux(iTracer,iCell,km1)
+ end do
+ end do
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+ do iTracer=1,grid % nTracers
+ s_update(iTracer,cell1,km0) = s_update(iTracer,cell1,km0) - &
+ h_flux(iTracer,iEdge) / grid % areaCell % array(cell1)
+ s_update(iTracer,cell2,km0) = s_update(iTracer,cell2,km0) + &
+ h_flux(iTracer,iEdge) / grid % areaCell % array(cell2)
+ end do
+ end if
+ end do
+
+ ! decouple from mass
+ if (k > 1) then
+ do iCell=1,grid % nCells
+ do iTracer=1,grid % nTracers
+ s_update(iTracer,iCell,km1) = s_update(iTracer,iCell,km1) / h_new(k-1,iCell)
+ end do
+ end do
+
+ do iCell=1,grid % nCells
+ do iTracer=1,grid % nTracers
+ scalar_new(iTracer,k-1,iCell) = s_update(iTracer,iCell,km1)
+ end do
+ end do
+ end if
+
+ ktmp = km1
+ km1 = km0
+ km0 = ktmp
+
+ end do
+
+ do iCell=1,grid % nCells
+ do iTracer=1,grid % nTracers
+ scalar_new(iTracer,grid % nVertLevels,iCell) = s_update(iTracer,iCell,km1) / h_new(grid%nVertLevels,iCell)
+ end do
+ end do
+
+ end subroutine advance_scalars_mono
+
+
subroutine compute_solve_diagnostics(dt, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
</font>
</pre>