<p><b>ringler@lanl.gov</b> 2011-03-30 16:16:30 -0600 (Wed, 30 Mar 2011)</p><p><br>
moved module_poisson_solver.F into the pvsw core<br>
<br>
tested the computation of poisson_weights<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/src/core_pvsw/Makefile
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Makefile 2011-03-24 22:25:27 UTC (rev 769)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Makefile 2011-03-30 22:16:30 UTC (rev 770)
@@ -3,6 +3,7 @@
OBJS = module_mpas_core.o \
module_test_cases.o \
module_advection.o \
+ module_poisson_solver.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -15,11 +16,13 @@
module_advection.o:
+module_poisson_solver.o:
+
module_time_integration.o:
module_global_diagnostics.o:
-module_mpas_core.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o module_advection.o
+module_mpas_core.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o module_advection.o module_poisson_solver.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-03-24 22:25:27 UTC (rev 769)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-03-30 22:16:30 UTC (rev 770)
@@ -14,6 +14,7 @@
use configure
use grid_types
use test_cases
+ use poisson_solver
implicit none
@@ -23,8 +24,19 @@
type (block_type), pointer :: block
+ write(6,*) 'calling sw_test_case'
if (.not. config_do_restart) call setup_sw_test_case(domain)
+ write(6,*) 'calling init_poisson'
+ call init_poisson(domain)
+
+ write(6,*) 'calling poisson'
+ call poisson(domain)
+
+ stop
+
+
+
!
! Initialize core
!
Added: branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F (rev 0)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-03-30 22:16:30 UTC (rev 770)
@@ -0,0 +1,242 @@
+module poisson_solver
+
+ use grid_types
+ use configure
+ use constants
+ use dmpar
+
+ implicit none
+
+ public :: init_poisson, gs_iteration, poisson
+
+ contains
+
+ subroutine init_poisson(domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: pre-compute coefficients used by poisson to invert the Laplacian
+ !
+ ! Input: grid meta data
+ !
+ ! Output: grid % poisson_weights - coefficients used in poisson solver
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ type (block_type), pointer :: block
+ integer :: j, iCell, jEdge, jCell
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ ! init weights to zero
+ block % mesh % poisson_weights % array(:,:) = 0.0
+
+ ! loop over all cells
+ do iCell=1,block % mesh % nCells
+ do j=1,block % mesh % nEdgesOnCell % array(iCell)
+
+ jEdge=block % mesh % edgesOnCell % array(j,iCell)
+ jCell=block % mesh % cellsOnCell % array(j,iCell)
+
+ block % mesh % poisson_weights % array(j,iCell) = block % mesh % dvEdge % array(jEdge) / &
+ block % mesh % dcEdge % array(jEdge)
+
+ write(10,*) iCell, j, jEdge, jCell, block % mesh % nEdgesOnCell % array(iCell)
+ write(10,*) block % mesh % poisson_weights % array(j,iCell)
+ write(10,*) block % mesh % dvEdge % array(jEdge)
+ write(10,*) block % mesh % dcEdge % array(jEdge)
+ write(10,*)
+
+ if(block % mesh % poisson_weights % array(j,iCell).le.0.0) write(6,*) ' too low '
+ if(block % mesh % poisson_weights % array(j,iCell).gt.1.0) write(6,*) ' too high'
+
+ enddo ! j
+ enddo ! iCell
+
+ block => block % next
+ enddo ! block
+
+
+ block => domain % blocklist
+ do while (associated(block))
+ write(6,*) ' min/max poisson_weights :', &
+ minval(block % mesh % poisson_weights % array(:,1:block % mesh % nCellsSolve)), &
+ maxval(block % mesh % poisson_weights % array(:,1:block % mesh % nCellsSolve))
+ block => block % next
+ enddo ! block
+
+ end subroutine init_poisson
+
+
+ subroutine poisson(domain)
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: iteration
+ logical :: converge
+ real (kind=RKIND) :: l2_psi_error, l2_chi_error
+ type (block_type), pointer :: block
+
+ iteration = 0
+ converge = .false.
+
+ do while (.not.converge)
+ iteration = iteration + 1
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ call gs_iteration(block % state % time_levs(2) % state, block % mesh)
+
+ if(mod(iteration,100).eq.0) call compute_poisson_error(block % state % time_levs(2) % state, block % mesh, l2_psi_error, l2_chi_error)
+ write(6,*) l2_psi_error, l2_chi_error
+
+! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % psi % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
+! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % chi % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
+ block => block % next
+
+ end do
+
+ end do
+
+ end subroutine poisson
+
+
+ subroutine gs_iteration(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: do one iteration using Guass Seidel to solve
+ ! psi and chi where del2(psi)=vorticity and del2(chi)=divergence
+ !
+ ! Input: niter: number of iterations to conduct
+ ! vorticity, divergence, initial guess for psi and chi
+ !
+ ! Output: update psi and chi
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ 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 :: vorticity, divergence, psi, chi
+ real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
+ real (kind=RKIND) :: sum1, sum2, sum3
+
+ cellsOnCell => grid % cellsOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ areaCell => grid % areaCell % array
+ nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
+
+ vorticity => state % vorticity % array
+ divergence => state % divergence % array
+ psi => state % psi % array
+ chi => state % vorticity % array
+
+ ! 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 + psi(k,jCell) * poisson_weights(j,iCell)
+ sum3 = sum3 + chi(k,jCell) * poisson_weights(j,iCell)
+ psi(k,iCell) = (vorticity(k,iCell)*areaCell(iCell) + sum2)/sum1
+ chi(k,iCell) = (divergence(k,iCell)*areaCell(iCell) + sum3)/sum1
+
+ enddo! loop over neighboring cells
+
+ enddo ! loop over cell centers
+ enddo ! loop over vertical levels
+
+ end subroutine gs_iteration
+
+
+
+ subroutine compute_poisson_error(state, grid, l2_psi_error, l2_chi_error)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: compute error in poisson solver
+ !
+ ! Input: psi, chi, vorticity and divergence
+ !
+ ! 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_psi_error, l2_chi_error
+
+ 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 :: vorticity, divergence, psi, chi
+ real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
+ real (kind=RKIND) :: sum_error_psi, sum_error_chi, sum_area, r1, r2
+
+ cellsOnCell => grid % cellsOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ areaCell => grid % areaCell % array
+ nCellsSolve = grid % nCellsSolve
+ nVertLevels = grid % nVertLevels
+
+ vorticity => state % vorticity % array
+ divergence => state % divergence % array
+ psi => state % psi % array
+ chi => state % vorticity % array
+
+
+ do k=1,nVertLevels
+
+ sum_error_psi = 0.0; sum_error_chi = 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) * (psi(k,jCell)-psi(k,iCell))
+ r2 = r2 + poisson_weights(j,iCell) * (chi(k,jCell)-chi(k,iCell))
+
+ enddo
+ r1 = r1 / areaCell(iCell)
+ r2 = r2 / areaCell(iCell)
+
+ sum_error_psi = sum_error_psi + areaCell(iCell)*(r1 - vorticity (k,iCell))**2
+ sum_error_chi = sum_error_chi + areaCell(iCell)*(r1 - divergence(k,iCell))**2
+ sum_area = sum_area + areaCell(iCell)
+
+ enddo
+
+ l2_psi_error = sqrt(sum_error_psi/sum_area)
+ l2_chi_error = sqrt(sum_error_chi/sum_area)
+
+ enddo
+
+ end subroutine compute_poisson_error
+
+end module poisson_solver
Modified: branches/pv_based_swm/mpas/src/operators/Makefile
===================================================================
--- branches/pv_based_swm/mpas/src/operators/Makefile 2011-03-24 22:25:27 UTC (rev 769)
+++ branches/pv_based_swm/mpas/src/operators/Makefile 2011-03-30 22:16:30 UTC (rev 770)
@@ -1,6 +1,6 @@
.SUFFIXES: .F .o
-OBJS = module_RBF_interpolation.o module_vector_reconstruction.o module_spline_interpolation.o module_poisson_solver.o
+OBJS = module_RBF_interpolation.o module_vector_reconstruction.o module_spline_interpolation.o
all: operators
@@ -10,7 +10,6 @@
module_vector_reconstruction.o: module_RBF_interpolation.o
module_RBF_interpolation.o:
module_spline_interpolation:
-module_poisson_solver.o:
clean:
$(RM) *.o *.mod *.f90 libops.a
Deleted: branches/pv_based_swm/mpas/src/operators/module_poisson_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/operators/module_poisson_solver.F 2011-03-24 22:25:27 UTC (rev 769)
+++ branches/pv_based_swm/mpas/src/operators/module_poisson_solver.F 2011-03-30 22:16:30 UTC (rev 770)
@@ -1,216 +0,0 @@
-module poisson_solver
-
- use grid_types
- use configure
- use constants
- use dmpar
-
- implicit none
-
- public :: init_poisson, gs_iteration, poisson
-
- contains
-
- subroutine init_poisson(grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: pre-compute coefficients used by poisson to invert the Laplacian
- !
- ! Input: grid meta data
- !
- ! Output: grid % poisson_weights - coefficients used in poisson solver
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
-
- integer :: nCells
- integer, dimension(:,:), pointer :: edgesOnCell
- integer, dimension(:), pointer :: nEdgesOnCell
- integer :: j, iCell, iEdge
- real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge
- real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
-
- !========================================================
- ! arrays filled and saved during init procedure
- !========================================================
- poisson_weights => grid % poisson_weights % array
- poisson_weights = 0.0
-
- !========================================================
- ! temporary variables needed for init procedure
- !========================================================
- dvEdge => grid % dvEdge % array
- dcEdge => grid % dcEdge % array
-
- ! loop over all cells
- do iCell=1,nCells
- do j=1,nEdgesOnCell(iCell)
- iEdge=edgesOnCell(j,iCell)
- poisson_weights(j,iCell) = dvEdge(iEdge)/dcEdge(iEdge)
- enddo ! j
- enddo ! iCell
-
- end subroutine init_poisson
-
-
- subroutine poisson(domain)
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: iteration
- logical :: converge
- type (block_type), pointer :: block
-
- iteration = 0
- converge = .false.
-
- do while (.not.converge)
- iteration = iteration + 1
-
- block => domain % blocklist
- do while (associated(block))
-
- call gs_iteration(1, block % state % time_levs(2) % state, block % mesh)
-
- if(mod(iteration,100).eq.0) call compute_poisson_error(block % state % time_levs(2) % state, block % mesh, converge)
-
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % psi % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % chi % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-
- block => block % next
-
- end do
-
- end do
-
- end subroutine poisson
-
-
- subroutine gs_iteration(nIter, state, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: do nIter iteration using Guass Seidel to solve
- ! psi and chi where del2(psi)=vorticity and del2(chi)=divergence
- !
- ! Input: niter: number of iterations to conduct
- ! vorticity, divergence, initial guess for psi and chi
- !
- ! Output: update psi and chi
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(inout) :: state
- type (mesh_type), intent(in) :: grid
-
- integer :: nIter, 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 :: vorticity, divergence, psi, chi
- real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
- real (kind=RKIND) :: sum1, sum2, sum3
-
- cellsOnCell => grid % cellsOnCell % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- areaCell => grid % areaCell % array
- nCells = grid % nCells
- nVertLevels = grid % nVertLevels
-
- vorticity => state % vorticity % array
- divergence => state % divergence % array
- psi => state % psi % array
- chi => state % vorticity % array
-
- ! do requested number of iterations
- do iter=1,nIter
-
- ! loop over the vertical
- do k = 1, nVertLevels
-
- ! loop over cell centers
- do iCell=1,nCells
-
- ! loop over neighboring cells
- do j=1,nEdgesOnCell(iCell)
- jCell = cellsOnCell(j,iCell)
-
- sum1 = sum1 + poisson_weights(j,iCell)
- sum2 = sum2 + psi(k,iCell) * poisson_weights(j,iCell)
- sum3 = sum3 + chi(k,iCell) * poisson_weights(j,iCell)
- psi(k,iCell) = (vorticity(k,iCell)*areaCell(iCell) + sum2)/sum1
- chi(k,iCell) = (divergence(k,iCell)*areaCell(iCell) + sum3)/sum1
-
- enddo! loop over neighboring cells
-
- enddo ! loop over cell centers
- enddo ! loop over vertical levels
- enddo ! loop over iterations
-
- end subroutine gs_iteration
-
-
-
- subroutine compute_poisson_error(state, grid, converge)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: compute error in poisson solver
- !
- ! Input: psi, chi, vorticity and divergence
- !
- ! Output: converge (true or false)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(in) :: state
- type (mesh_type), intent(in) :: grid
- logical, intent(out) :: converge
-
- 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 :: vorticity, divergence, psi, chi
- real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
- real (kind=RKIND) :: error_psi, error_chi
-
- cellsOnCell => grid % cellsOnCell % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- areaCell => grid % areaCell % array
- nCellsSolve = grid % nCellsSolve
- nVertLevels = grid % nVertLevels
-
- vorticity => state % vorticity % array
- divergence => state % divergence % array
- psi => state % psi % array
- chi => state % vorticity % array
-
-
- error_psi = 0.0; error_chi = 0.0
- do iCell=1,nCellsSolve
-
-
-
-
-
-
- enddo
-
-
-
-
- converge = .false.
-
-
- end subroutine compute_poisson_error
-
-end module poisson_solver
</font>
</pre>