<p><b>ringler@lanl.gov</b> 2011-04-21 13:54:24 -0600 (Thu, 21 Apr 2011)</p><p><br>
add module for alpha closure<br>
</p><hr noshade><pre><font color="gray">Added: branches/pv_based_swm/mpas/src/core_pvsw/module_alpha.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_alpha.F (rev 0)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_alpha.F 2011-04-21 19:54:24 UTC (rev 799)
@@ -0,0 +1,273 @@
+module alpha_solver
+
+ use grid_types
+ use configure
+ use constants
+ use poisson_solver
+ use dmpar
+
+ implicit none
+
+ real (KIND=RKIND), parameter :: threshold = 1.0e-8
+ integer, parameter :: nIterPerGS = 2
+
+ public :: init_alpha, gs_iteration_alpha, alpha
+
+ contains
+
+ subroutine init_alpha(domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: put any init stuff in here
+ !
+ ! Input: grid meta data
+ !
+ ! Output:
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ type (block_type), pointer :: block
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ block => block % next
+ enddo ! block
+
+ end subroutine init_alpha
+
+
+ subroutine alpha(domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: given the rough vorticity and divergence, find the smooth vorticity and divergence
+ !
+ ! Input: grid meta data
+ !
+ ! Output: vorSmooth, divSmooth
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: iteration
+ logical :: converge
+ real (kind=RKIND) :: l2_vorSmooth_error, l2_divSmooth_error, globalVorSmoothError, globalDivSmoothError, globalArea, localArea, rflag, gflag
+ type (block_type), pointer :: block
+
+ iteration = 0; l2_vorSmooth_error = 1.0; l2_divSmooth_error = 1.0
+ converge = .false.
+
+ ! iterate until we converge
+ do while (.not.converge)
+ iteration = iteration + 1
+
+ ! do a few Guass-Seidel iterations
+ block => domain % blocklist
+ do while (associated(block))
+ call gs_iteration_alpha(block % state % time_levs(2) % state, block % mesh)
+ block => block % next
+ end do
+
+ ! do a boundary update on vorSmooth and divSmooth
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorSmooth % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ if(.not.config_bve) then
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divSmooth % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ endif
+ block => block % next
+ end do
+
+
+ ! check the global error every so often
+ if(mod(iteration,100).eq.0) then
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_alpha_error(block % state % time_levs(2) % state, block % mesh, l2_vorSmooth_error, l2_divSmooth_error, localArea)
+ call dmpar_sum_real(domain % dminfo, l2_vorSmooth_error, globalVorSmoothError)
+ call dmpar_sum_real(domain % dminfo, l2_divSmooth_error, globalDivSmoothError)
+ call dmpar_sum_real(domain % dminfo, localArea, globalArea)
+ block => block % next
+ end do
+ l2_vorSmooth_error = sqrt(globalVorSmoothError/globalArea)
+ l2_divSmooth_error = sqrt(globalDivSmoothError/globalArea)
+ if(l2_vorSmooth_error.lt.threshold.and.l2_divSmooth_error.lt.threshold) converge=.true.
+ if(iteration.gt.10000) converge=.true.
+ endif
+
+
+ ! make sure that everyone agrees on convergence
+ if(mod(iteration,100).eq.0) then
+ block => domain % blocklist
+ do while (associated(block))
+ rflag = 0.0
+ if(converge) rflag=1.0
+ call dmpar_min_real(domain % dminfo, rflag, gflag)
+ block => block % next
+ end do
+ converge=.false.
+ if(gflag.gt.0.5) converge=.true.
+ write(6,10) iteration, l2_vorSmooth_error, l2_divSmooth_error, rflag, gflag
+ 10 format(i8,6e20.10)
+ endif
+
+ enddo ! do while not converged
+
+ ! copy vorSmooth and divSmooth into the other time level
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % vorSmooth % array(:,:) = block % state % time_levs(2) % state % vorSmooth % array(:,:)
+ block % state % time_levs(1) % state % divSmooth % array(:,:) = block % state % time_levs(2) % state % divSmooth % array(:,:)
+ block => block % next
+ end do
+
+ end subroutine alpha
+
+
+ subroutine gs_iteration_alpha(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: do one or more iteration(s) using Guass Seidel to solve
+ ! vorSmooth and divSmooth where del2(vorSmooth)=vor and del2(divSmooth)=div
+ !
+ ! Input: vor, div, initial guess for vorSmooth and divSmooth
+ !
+ ! Output: update vorSmooth and divSmooth
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ integer :: nCells
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnCell
+ integer :: iCell, jCell, iter, k, j, nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: vor, div, vorSmooth, divSmooth
+ real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
+ real (kind=RKIND) :: sum1, sum2, sum3
+ real (kind=RKIND), parameter :: sor = 1.5
+
+ cellsOnCell => grid % cellsOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ areaCell => grid % areaCell % array
+ poisson_weights => grid % poisson_weights % array
+ nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
+
+ vor => state % vor % array
+ div => state % div % array
+ vorSmooth => state % vorSmooth % array
+ divSmooth => state % divSmooth % array
+
+ ! do nIter GS iterations
+ do iter=1,nIterPerGS
+
+ ! loop over the vertical
+ do k = 1, nVertLevels
+
+ ! loop over cell centers
+ do iCell=1,nCells
+ sum1=0.0; sum2=0.0; sum3=0.0
+
+ ! loop over neighboring cells
+ do j=1,nEdgesOnCell(iCell)
+ jCell = cellsOnCell(j,iCell)
+
+ sum1 = sum1 + poisson_weights(j,iCell)
+ sum2 = sum2 + poisson_weights(j,iCell) * vorSmooth(k,jCell)
+ sum3 = sum3 + poisson_weights(j,iCell) * divSmooth(k,jCell)
+
+ enddo ! loop over edges
+
+ sum1 = sum1 * alphasq / areaCell(iCell)
+ sum2 = sum2 * alphasq / areaCell(iCell)
+ sum3 = sum3 * alphasq / areaCell(iCell)
+
+ vorSmooth(k,iCell) = (1.0-sor)*vorSmooth(k,iCell) + sor*(sum2 + vor(k,iCell))/(1.0+sum1)
+ divSmooth(k,iCell) = (1.0-sor)*divSmooth(k,iCell) + sor*(sum3 + div(k,iCell))/(1.0+sum1)
+
+ enddo ! loop over cell centers
+ enddo ! loop over vertical levels
+
+ enddo ! loop over iterations
+
+ end subroutine gs_iteration_alpha
+
+
+ subroutine compute_alpha_error(state, grid, l2_vorSmooth_error, l2_divSmooth_error, localArea)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: compute error in poisson solver
+ !
+ ! Input: vorSmooth, divSmooth, vor and div
+ !
+ ! Output: converge (true or false)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(in) :: state
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), intent(out) :: l2_vorSmooth_error, l2_divSmooth_error, localArea
+
+ integer :: nCellsSolve
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnCell
+ integer :: iCell, jCell, iter, k, j, nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: vor, div, vorSmooth, divSmooth
+ real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
+ real (kind=RKIND) :: sum_error_vorSmooth, sum_error_divSmooth, sum_area, r1, r2
+
+ poisson_weights => grid % poisson_weights % array
+ cellsOnCell => grid % cellsOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ areaCell => grid % areaCell % array
+ nCellsSolve = grid % nCellsSolve
+ nVertLevels = grid % nVertLevels
+
+ vor => state % vor % array
+ div => state % div % array
+ vorSmooth => state % vorSmooth % array
+ divSmooth => state % divSmooth % array
+
+ do k=1,nVertLevels
+
+ sum_error_vorSmooth = 0.0; sum_error_divSmooth = 0.0; sum_area = 0.0
+ do iCell=1,nCellsSolve
+
+ r1=0.0; r2=0.0
+ do j=1,nEdgesOnCell(iCell)
+
+ jCell = cellsOnCell(j,iCell)
+ r1 = r1 + poisson_weights(j,iCell) * (vorSmooth(k,jCell)-vorSmooth(k,iCell))
+ r2 = r2 + poisson_weights(j,iCell) * (divSmooth(k,jCell)-divSmooth(k,iCell))
+
+ enddo
+ r1 = vorSmooth(k,iCell) - r1 * alphasq / areaCell(iCell)
+ r2 = divSmooth(k,iCell) - r2 * alphasq / areaCell(iCell)
+
+ sum_error_vorSmooth = sum_error_vorSmooth + areaCell(iCell)*(r1 - vor (k,iCell))**2
+ sum_error_divSmooth = sum_error_divSmooth + areaCell(iCell)*(r2 - div(k,iCell))**2
+ sum_area = sum_area + areaCell(iCell)
+
+ enddo
+
+ l2_vorSmooth_error = sum_error_vorSmooth
+ l2_divSmooth_error = sum_error_divSmooth
+ localArea = sum_area
+
+ enddo
+
+ end subroutine compute_alpha_error
+
+end module alpha_solver
</font>
</pre>