<p><b>ringler@lanl.gov</b> 2011-03-30 20:35:45 -0600 (Wed, 30 Mar 2011)</p><p><br>
poisson solver working<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-03-30 22:16:30 UTC (rev 770)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-03-31 02:35:45 UTC (rev 771)
@@ -25,7 +25,8 @@
type (domain_type), intent(inout) :: domain
type (block_type), pointer :: block
- integer :: j, iCell, jEdge, jCell
+ integer :: k, j, iCell, jEdge, jCell, kEdge, kCell
+ real (KIND=RKIND) :: r,s
block => domain % blocklist
do while (associated(block))
@@ -42,7 +43,8 @@
block % mesh % poisson_weights % array(j,iCell) = block % mesh % dvEdge % array(jEdge) / &
block % mesh % dcEdge % array(jEdge)
-
+
+ ! dump weights
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)
@@ -59,6 +61,7 @@
enddo ! block
+ ! write out min and max
block => domain % blocklist
do while (associated(block))
write(6,*) ' min/max poisson_weights :', &
@@ -67,6 +70,42 @@
block => block % next
enddo ! block
+ ! check for symmetry
+ block => domain % blocklist
+ do while (associated(block))
+ do iCell=1,block % mesh % nCellsSolve
+ do j=1,block % mesh % nEdgesOnCell % array(iCell)
+ r = block % mesh % poisson_weights % array(j,iCell)
+
+ jEdge=block % mesh % edgesOnCell % array(j,iCell)
+ jCell=block % mesh % cellsOnCell % array(j,iCell)
+
+ do k=1,block % mesh % nEdgesOnCell % array(jCell)
+ s = block % mesh % poisson_weights % array(k,jCell)
+
+ kEdge=block % mesh % edgesOnCell % array(k,jCell)
+ kCell=block % mesh % cellsOnCell % array(k,jCell)
+
+ if(kEdge.eq.jEdge) then
+ if(kCell.ne.iCell) then
+ write(6,*) ' kCell and iCell do not match'
+ write(6,*) iCell, kCell
+ endif
+ if(abs(s-r).gt.1.0e-10) then
+ write(6,*) ' weights not symmetric'
+ write(6,*) r,s
+ endif
+ endif
+
+ enddo
+
+ enddo
+ enddo
+
+ block => block % next
+ enddo ! block
+
+
end subroutine init_poisson
@@ -84,6 +123,32 @@
iteration = 0
converge = .false.
+ write(6,*) ' starting poisson '
+ write(6,*) ' init psi,chi,vorticity,divergence'
+ block => domain % blocklist
+ do while (associated(block))
+
+ block % state % time_levs(2) % state % psi % array(:,:) = 0.0
+ block % state % time_levs(2) % state % chi % array(:,:) = 0.0
+ block % state % time_levs(2) % state % vorticity % array(:,:) = 0.0
+ block % state % time_levs(2) % state % divergence % array(:,:) = 0.0
+
+ block => block % next
+
+ end do
+
+ write(6,*) ' set vorticity'
+ block => domain % blocklist
+ do while (associated(block))
+
+ block % state % time_levs(2) % state % vorticity % array(1,1225) = 1.0e-4
+ block % state % time_levs(2) % state % vorticity % array(:,:) = block % state % time_levs(2) % state % vorticity % array(:,:) - 1.0e-4/2500
+
+ block => block % next
+
+ end do
+
+ write(6,*) 'start iteration'
do while (.not.converge)
iteration = iteration + 1
@@ -93,7 +158,7 @@
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
+ if(mod(iteration,100).eq.0) 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, &
@@ -117,8 +182,7 @@
! 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
+ ! Input: vorticity, divergence, initial guess for psi and chi
!
! Output: update psi and chi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -128,7 +192,7 @@
type (state_type), intent(inout) :: state
type (mesh_type), intent(in) :: grid
- integer :: nCells
+ integer :: nCellsSolve
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: cellsOnCell
integer :: iCell, jCell, iter, k, j, nVertLevels
@@ -137,11 +201,12 @@
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
+ cellsOnCell => grid % cellsOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ areaCell => grid % areaCell % array
+ nCellsSolve = grid % nCellsSolve
+ nVertLevels = grid % nVertLevels
+ poisson_weights => grid % poisson_weights % array
vorticity => state % vorticity % array
divergence => state % divergence % array
@@ -152,7 +217,7 @@
do k = 1, nVertLevels
! loop over cell centers
- do iCell=1,nCells
+ do iCell=1,nCellsSolve
sum1=0.0; sum2=0.0; sum3=0.0
! loop over neighboring cells
@@ -160,13 +225,14 @@
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
+ sum2 = sum2 + poisson_weights(j,iCell) * psi(k,jCell)
+ sum3 = sum3 + poisson_weights(j,iCell) * chi(k,jCell)
- enddo! loop over neighboring cells
+ enddo ! loop over edges
+ psi(k,iCell) = (sum2 - vorticity (k,iCell)*areaCell(iCell))/sum1
+ chi(k,iCell) = (sum3 - divergence(k,iCell)*areaCell(iCell))/sum1
+
enddo ! loop over cell centers
enddo ! loop over vertical levels
@@ -198,11 +264,12 @@
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
+ poisson_weights => grid % poisson_weights % array
+ 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
@@ -227,7 +294,7 @@
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_error_chi = sum_error_chi + areaCell(iCell)*(r2 - divergence(k,iCell))**2
sum_area = sum_area + areaCell(iCell)
enddo
</font>
</pre>