<p><b>ringler@lanl.gov</b> 2011-03-31 16:00:05 -0600 (Thu, 31 Mar 2011)</p><p><br>
added poisson solver to the time integration scheme<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-31 02:35:45 UTC (rev 771)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Makefile 2011-03-31 22:00:05 UTC (rev 772)
@@ -18,7 +18,7 @@
module_poisson_solver.o:
-module_time_integration.o:
+module_time_integration.o: module_poisson_solver.o
module_global_diagnostics.o:
Modified: branches/pv_based_swm/mpas/src/core_pvsw/Registry
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-03-31 02:35:45 UTC (rev 771)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-03-31 22:00:05 UTC (rev 772)
@@ -122,8 +122,9 @@
var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
var persistent real h ( nVertLevels nCells Time ) 2 iro h state - -
var persistent real tracers ( nTracers nVertLevels nCells Time ) 2 iro tracers state - -
-var persistent real pv ( nVertLevels nCells Time ) vorticity iro pv state - -
-var persistent real div ( nVertLevels nCells Time ) divergence iro div state - -
+var persistent real pv ( nVertLevels nCells Time ) 2 iro pv state - -
+var persistent real vor ( nVertLevels nCells Time ) 2 iro vor state - -
+var persistent real div ( nVertLevels nCells Time ) 2 iro div state - -
var persistent real psi ( nVertLevels nCells Time ) 2 iro psi state - -
var persistent real chi ( nVertLevels nCells Time ) 2 iro chi state - -
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-31 02:35:45 UTC (rev 771)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-03-31 22:00:05 UTC (rev 772)
@@ -33,10 +33,6 @@
write(6,*) 'calling poisson'
call poisson(domain)
- stop
-
-
-
!
! Initialize core
!
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-31 02:35:45 UTC (rev 771)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-03-31 22:00:05 UTC (rev 772)
@@ -7,6 +7,8 @@
implicit none
+ real (KIND=RKIND), parameter :: threshold = 1.0e-7
+
public :: init_poisson, gs_iteration, poisson
contains
@@ -120,69 +122,64 @@
real (kind=RKIND) :: l2_psi_error, l2_chi_error
type (block_type), pointer :: block
- iteration = 0
+ iteration = 0; l2_psi_error = 1.0; l2_chi_error = 1.0
converge = .false.
- write(6,*) ' starting poisson '
- write(6,*) ' init psi,chi,vorticity,divergence'
- block => domain % blocklist
- do while (associated(block))
+! update vor and div to get started
+! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vor % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- 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
+! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % div % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- 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
+ iteration = iteration + 1
- block => domain % blocklist
- do while (associated(block))
+ block => domain % blocklist
+ do while (associated(block))
- call gs_iteration(block % state % time_levs(2) % state, block % mesh)
+ 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)
- if(mod(iteration,100).eq.0) call compute_poisson_error(block % state % time_levs(2) % state, block % mesh, 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, &
+! 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)
-! 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)
+! sum error
-! 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
- block => block % next
+ end do
- end do
+ if(l2_psi_error.lt.threshold.and.l2_chi_error.lt.threshold) converge=.true.
+ if(iteration.gt.10000) converge=.true.
+ block => domain % blocklist
+ do while (associated(block))
+ if(mod(iteration,100).eq.0) write(6,10) iteration, l2_psi_error, l2_chi_error, maxval(block % state % time_levs(2) % state % chi % array(:,:))
+ block => block % next
+ end do
+
+ 10 format(i8,3e20.10)
end do
+ write(6,*) ' exiting poisson'
+
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
+ ! psi and chi where del2(psi)=vor and del2(chi)=div
!
- ! Input: vorticity, divergence, initial guess for psi and chi
+ ! Input: vor, div, initial guess for psi and chi
!
! Output: update psi and chi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -197,7 +194,7 @@
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 :: vor, div, psi, chi
real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
real (kind=RKIND) :: sum1, sum2, sum3
@@ -208,10 +205,10 @@
nVertLevels = grid % nVertLevels
poisson_weights => grid % poisson_weights % array
- vorticity => state % vorticity % array
- divergence => state % divergence % array
+ vor => state % vor % array
+ div => state % div % array
psi => state % psi % array
- chi => state % vorticity % array
+ chi => state % vor % array
! loop over the vertical
do k = 1, nVertLevels
@@ -230,8 +227,8 @@
enddo ! loop over edges
- psi(k,iCell) = (sum2 - vorticity (k,iCell)*areaCell(iCell))/sum1
- chi(k,iCell) = (sum3 - divergence(k,iCell)*areaCell(iCell))/sum1
+ psi(k,iCell) = (sum2 - vor (k,iCell)*areaCell(iCell))/sum1
+ chi(k,iCell) = (sum3 - div(k,iCell)*areaCell(iCell))/sum1
enddo ! loop over cell centers
enddo ! loop over vertical levels
@@ -239,12 +236,11 @@
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
+ ! Input: psi, chi, vor and div
!
! Output: converge (true or false)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -260,7 +256,7 @@
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 :: vor, div, psi, chi
real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
real (kind=RKIND) :: sum_error_psi, sum_error_chi, sum_area, r1, r2
@@ -271,10 +267,10 @@
nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
- vorticity => state % vorticity % array
- divergence => state % divergence % array
+ vor => state % vor % array
+ div => state % div % array
psi => state % psi % array
- chi => state % vorticity % array
+ chi => state % vor % array
do k=1,nVertLevels
@@ -293,8 +289,8 @@
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)*(r2 - divergence(k,iCell))**2
+ sum_error_psi = sum_error_psi + areaCell(iCell)*(r1 - vor (k,iCell))**2
+ sum_error_chi = sum_error_chi + areaCell(iCell)*(r2 - div(k,iCell))**2
sum_area = sum_area + areaCell(iCell)
enddo
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_test_cases.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_test_cases.F 2011-03-31 02:35:45 UTC (rev 771)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_test_cases.F 2011-03-31 22:00:05 UTC (rev 772)
@@ -242,7 +242,14 @@
sin(grid%latVertex%array(iVtx)) * cos(alpha) &
)
end do
+ do iCell=1,grid % nCells
+ grid % fCell % array(iCell) = 2.0 * omega * &
+ (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
+ sin(grid%latCell%array(iCell)) * cos(alpha) &
+ )
+ end do
+
!
! Initialize height field (actually, fluid thickness field)
!
@@ -255,6 +262,22 @@
gravity
end do
+ ! Initialize psi, chi, pv, vor and div (pv solver)
+ do iCell=1,grid % nCells
+ state % psi % array(1,iCell) = -a * u0 * ( &
+ sin(grid%latCell%array(iCell)) * cos(alpha) - &
+ cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) &
+ )
+ state % chi % array(1,iCell) = 0.0
+ state % div % array(1,iCell) = 0.0
+ state % pv % array(1,iCell) = ( (2.0*u0/a + 2.0*omega)* &
+ (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
+ sin(grid%latCell%array(iCell)) * cos(alpha) ) + grid%fCell%array(iCell) ) / &
+ state % h % array(1,iCell)
+ state % vor % array(1,iCell) = state % h % array(1,iCell) * state % pv % array(1,iCell) - grid%fCell%array(iCell)
+ end do
+
+
end subroutine sw_test_case_2
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-03-31 02:35:45 UTC (rev 771)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-03-31 22:00:05 UTC (rev 772)
@@ -4,6 +4,7 @@
use grid_types
use configure
use constants
+ use poisson_solver
use dmpar
@@ -159,6 +160,16 @@
block => block % next
end do
+! --- invert elliptic equations for psi and chi
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(2) % state % vor % array(:,:) = provis % vorticity_cell % array(:,:)
+ block % state % time_levs(2) % state % div % array(:,:) = provis % divergence % array(:,:)
+ block => block % next
+ end do
+ write(11,*) 'calling poisson from module_time_integration'
+ call poisson(domain)
+
! --- compute next substep state
if (rk_step < 4) then
</font>
</pre>