<p><b>qchen3@fsu.edu</b> 2012-01-23 11:10:00 -0700 (Mon, 23 Jan 2012)</p><p>BRANCH COMMIT<br>
<br>
To prepare for developing a new test case with a double-periodic domain and propogating gravity waves.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/fvswe/Makefile
===================================================================
--- branches/fvswe/Makefile 2012-01-23 15:54:21 UTC (rev 1408)
+++ branches/fvswe/Makefile 2012-01-23 18:10:00 UTC (rev 1409)
@@ -56,6 +56,18 @@
"CORE = $(CORE)" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+pgi-cray:
+ ( make all \
+ "FC = ftn" \
+ "CC = cc" \
+ "SFC = ftn" \
+ "SCC = cc" \
+ "FFLAGS = -r8 -O3 -byteswapio" \
+ "CFLAGS = -O3" \
+ "LDFLAGS = -O3" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+
pgi-llnl:
( make all \
"FC = mpipgf90" \
Modified: branches/fvswe/namelist.input
===================================================================
--- branches/fvswe/namelist.input 2012-01-23 15:54:21 UTC (rev 1408)
+++ branches/fvswe/namelist.input 2012-01-23 18:10:00 UTC (rev 1409)
@@ -1 +1 @@
-link namelist.input.sw
\ No newline at end of file
+link namelist.input.swe
\ No newline at end of file
Modified: branches/fvswe/namelist.input.sw
===================================================================
--- branches/fvswe/namelist.input.sw 2012-01-23 15:54:21 UTC (rev 1408)
+++ branches/fvswe/namelist.input.sw 2012-01-23 18:10:00 UTC (rev 1409)
@@ -2,14 +2,15 @@
config_test_case = 5
config_time_integration = 'RK4'
config_dt = 172.8
- config_ntimesteps = 7500
- config_output_interval = 500
+ config_ntimesteps = 25000
+ config_output_interval = 2500
config_stats_interval = 0
config_h_ScaleWithMesh = .false.
config_h_mom_eddy_visc2 = 0.0
- config_h_mom_eddy_visc4 = 0.0
+ config_h_mom_eddy_visc4 = 10000000000000.0
config_h_tracer_eddy_diff2 = 0.0
config_h_tracer_eddy_diff4 = 0.0
+ config_apvm_upwinding = 0.0
config_thickness_adv_order = 2
config_tracer_adv_order = 2
config_positive_definite = .false.
Modified: branches/fvswe/namelist.input.swe
===================================================================
--- branches/fvswe/namelist.input.swe 2012-01-23 15:54:21 UTC (rev 1408)
+++ branches/fvswe/namelist.input.swe 2012-01-23 18:10:00 UTC (rev 1409)
@@ -2,12 +2,13 @@
config_test_case = 5
config_time_integration = 'RK4'
config_dt = 172.8
- config_ntimesteps = 7500
- config_output_interval = 500
+ config_ntimesteps = 25000
+ config_output_interval = 2500
config_stats_interval = 0
config_h_ScaleWithMesh = .false.
config_h_mom_eddy_visc2 = 0.0
config_h_mom_eddy_visc4 = 0.0
+ config_apvm_upwinding = 0.0020
config_h_tracer_eddy_diff2 = 0.0
config_h_tracer_eddy_diff4 = 0.0
config_thickness_adv_order = 2
Modified: branches/fvswe/src/core_swe/module_test_cases.F
===================================================================
--- branches/fvswe/src/core_swe/module_test_cases.F 2012-01-23 15:54:21 UTC (rev 1408)
+++ branches/fvswe/src/core_swe/module_test_cases.F 2012-01-23 18:10:00 UTC (rev 1409)
@@ -213,7 +213,10 @@
! Initialize wind field
!
allocate(psiVertex(grid % nVertices))
- allocate(psiCell(grid % nCells))
+ allocate(psiCell(grid % nCells+1))
+ ! psiCell(nCells+1) is a dummy cell
+ psiCell(grid % nCells + 1) = 0.0
+
do iVtx=1,grid % nVertices
psiVertex(iVtx) = -a * u0 * ( &
sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
@@ -335,7 +338,10 @@
! Initialize wind field
!
allocate(psiVertex(grid % nVertices))
- allocate(psiCell(grid % nCells))
+ allocate(psiCell(grid % nCells + 1))
+ ! psiCell(nCells+1) is a dummy cell
+ psiCell(grid % nCells + 1) = 0.0
+
do iVtx=1,grid % nVertices
psiVertex(iVtx) = -a * u0 * ( &
sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
@@ -349,7 +355,6 @@
)
end do
-
do iEdge=1,grid % nEdges
state % u % array(1,iEdge) = -1.0 * ( &
psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
Modified: branches/fvswe/src/core_swe/module_time_integration.F
===================================================================
--- branches/fvswe/src/core_swe/module_time_integration.F 2012-01-23 15:54:21 UTC (rev 1408)
+++ branches/fvswe/src/core_swe/module_time_integration.F 2012-01-23 18:10:00 UTC (rev 1409)
@@ -82,7 +82,6 @@
!
block => domain % blocklist
do while (associated(block))
-
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
block % state % time_levs(2) % state % v % array(:,:) = block % state % time_levs(1) % state % v % array(:,:)
block % state % time_levs(2) % state % h_cell % array(:,:) = block % state % time_levs(1) % state % h_cell % array(:,:)
@@ -125,18 +124,36 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, provis % h_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+! call dmpar_exch_halo_field2dReal(domain % dminfo, provis % u % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nEdges, &
+! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+! call dmpar_exch_halo_field2dReal(domain % dminfo, provis % v % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nEdges, &
+! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % ke_vertex % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+! call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_vertex % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nVertices, &
+! block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % ke_cell % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_cell % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
if (config_h_mom_eddy_visc4 > 0.0) then
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence_cell % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % divergence_cell % 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 % divergence_vertex % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % divergence_vertex % array(:,:), &
block % mesh % nVertLevels, block % mesh % nVertices, &
block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity_vertex % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % vorticity_vertex % array(:,:), &
block % mesh % nVertLevels, block % mesh % nVertices, &
block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity_cell % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % vorticity_cell % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
end if
@@ -194,8 +211,8 @@
provis % tracers % array(:,k,iCell) = ( &
block % state % time_levs(1) % state % h_cell % array(k,iCell) * &
block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
- ) / provis % h_cell % array(k,iCell)
+ + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
+ ) / provis % h_cell % array(k,iCell)
end do
end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
@@ -302,9 +319,7 @@
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
- !real (kind=RKIND) :: ke_edge
-
h_cell => s % h_cell % array
h_vertex => s % h_vertex % array
h_edge => s % h_edge % array
@@ -316,12 +331,10 @@
vorticity_vertex => s % vorticity_vertex % array
divergence_cell => s % divergence_cell % array
divergence_vertex => s % divergence_vertex % array
- !ke_edge => s % ke_edge % array
ke_cell => s % ke_cell % array
ke_vertex => s % ke_vertex % array
pv_edge => s % pv_edge % array
- !weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
cellsOnEdge => grid % cellsOnEdge % array
cellsOnVertex => grid % cellsOnVertex % array
@@ -362,7 +375,7 @@
!
tend_h_cell(:,:) = 0.0
tend_h_vertex(:,:) = 0.0
- do iEdge=1,nEdgesSolve
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
@@ -380,13 +393,13 @@
end do
end do
- do iCell=1,grid % nCellsSolve
+ do iCell=1,grid % nCells
do k=1,nVertLevels
tend_h_cell(k,iCell) = tend_h_cell(k,iCell) / areaCell(iCell)
end do
end do
- do iVertex=1,grid % nVerticesSolve
+ do iVertex=1,grid % nVertices
do k=1,nVertLevels
tend_h_vertex(k,iVertex) = tend_h_vertex(k,iVertex) / areaTriangle(iVertex)
end do
@@ -398,7 +411,7 @@
tend_u(:,:) = 0.0
tend_v(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
@@ -430,12 +443,12 @@
do k=1,nVertLevels
u_diffusion = ( divergence_cell(k,cell2) - divergence_cell(k,cell1) ) / dcEdge(iEdge) &
-(vorticity_vertex(k,vertex2) - vorticity_vertex(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
+ u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
v_diffusion = ( divergence_vertex(k,vertex2) - divergence_vertex(k,vertex1) ) / dvEdge(iEdge) &
+ (vorticity_cell(k,cell2) - vorticity_cell(k,cell1) ) / dcEdge(iEdge)
- v_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * v_diffusion
+ v_diffusion = config_h_mom_eddy_visc2 * v_diffusion
tend_v(k,iEdge) = tend_v(k,iEdge) + v_diffusion
end do
end do
@@ -568,8 +581,8 @@
+ ( delsq_vorticity_cell(k,cell2) &
- delsq_vorticity_cell(k,cell1) ) / dcEdge(iEdge)
- u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
- v_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * v_diffusion
+ u_diffusion = config_h_mom_eddy_visc4 * u_diffusion
+ v_diffusion = config_h_mom_eddy_visc4 * v_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
tend_v(k,iEdge) = tend_v(k,iEdge) - v_diffusion
@@ -636,7 +649,7 @@
gradPVn, gradPVt, divergence_cell, divergence_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r, h1, h2
+ real (kind=RKIND) :: r, h1, h2, m_ke, m_h, m_Jacobian, Jacobian, coef, alpha
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
@@ -691,6 +704,8 @@
boundaryEdge => grid % boundaryEdge % array
boundaryCell => grid % boundaryCell % array
+ alpha = config_apvm_upwinding
+
!
! Find those cells that have an edge on the boundary
!
@@ -729,6 +744,7 @@
! used to when reading for edges that do not exist
!
u(:,nEdges+1) = 0.0
+ v(:,nEdges+1) = 0.0
!
! Compute circulation and relative vorticity at each vertex
@@ -736,7 +752,7 @@
circulation_vertex(:,:) = 0.0
circulation_cell(:,:) = 0.0
- do iEdge=1,nEdgesSolve
+ do iEdge=1,nEdges
do k=1,nVertLevels
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
@@ -750,13 +766,13 @@
end do
end do
- do iVertex=1,nVerticesSolve
+ do iVertex=1,nVertices
do k=1,nVertLevels
vorticity_vertex(k,iVertex) = circulation_vertex(k,iVertex) / areaTriangle(iVertex)
end do
end do
- do iCell = 1, nCellsSolve
+ do iCell = 1, nCells
do k = 1, nVertLevels
vorticity_cell(k,iCell) = circulation_cell(k,iCell) / areaCell(iCell)
end do
@@ -851,19 +867,19 @@
! Compute pv at vertices and cells, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
- do iVertex = 1,nVerticesSolve
+ do iVertex = 1,nVertices
do k=1,nVertLevels
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity_vertex(k,iVertex)) / h_vertex(k,iVertex)
end do
end do
- do iCell = 1, nCellsSolve
+ do iCell = 1, nCells
do k = 1, nVertLevels
pv_cell(k,iCell) = (fCell(iCell) + vorticity_cell(k,iCell)) / h_cell(k,iCell)
enddo
end do
- do iEdge = 1, nEdgesSolve
+ do iEdge = 1, nEdges
do k = 1, nVertLevels
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -878,6 +894,7 @@
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
+ gradPVt(:,:) = 0.0
do iEdge = 1,nEdges
do k = 1,nVertLevels
gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
@@ -886,16 +903,6 @@
enddo
!
- ! Modify PV edge with upstream bias.
- !
-! do iEdge = 1,nEdges
-! do k = 1,nVertLevels
-! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
-! enddo
-! enddo
-
-
- !
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
@@ -909,15 +916,65 @@
endif
enddo
+ !
! Modify PV edge with upstream bias.
!
-! do iEdge = 1,nEdges
-! do k = 1,nVertLevels
-! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
-! enddo
-! enddo
+ !do iEdge = 1,nEdges
+ ! do k = 1,nVertLevels
+ ! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ ! enddo
+ !enddo
+ !
+ ! Modify PV edge with upstream bias.
+ !
+ !do iEdge = 1,nEdges
+ ! do k = 1,nVertLevels
+ ! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ ! enddo
+ !enddo
+ !
+ ! Scale-aware APVM
+ !
+ if (config_apvm_upwinding > epsilon(1D0)) then
+ do k = 1, nVertLevels
+ m_ke = sum(ke_cell(k,1:nCells) * areaCell(1:nCells)) / sum(areaCell(1:nCells))
+ m_ke = m_ke + sum(ke_vertex(k,1:nVertices) * areaTriangle(1:nVertices)) / sum(areaTriangle(1:nVertices))
+ m_ke = 0.5 * m_ke
+
+ m_h = sum(h_cell(k,1:nCells) * areaCell(1:nCells)) / sum(areaCell(1:nCells))
+ m_h = m_h + sum(h_vertex(k,1:nVertices) * areaTriangle(1:nVertices)) / sum(areaTriangle(1:nVertices))
+ m_h = 0.5 * m_h
+
+ m_Jacobian = 0.0
+ do iEdge = 1, nEdges
+ !m_Jacobian = m_Jacobian + abs(u(k,iEdge)*gradPVn(k,iEdge) + v(k,iEdge)*gradPVt(k,iEdge))
+ m_Jacobian = m_Jacobian + (u(k,iEdge)*gradPVn(k,iEdge) + v(k,iEdge)*gradPVt(k,iEdge))**2 * dcEdge(iEdge) * dvEdge(iEdge)
+ end do
+ m_Jacobian = 0.5 * m_Jacobian / sum(areaCell(1:nCells))
+ m_Jacobian = sqrt(m_Jacobian)
+! write(*,*) 'm_h = ', m_h
+! write(*,*) 'm_ke = ', m_ke
+! write(*,*) 'm_Jacobian = ', m_Jacobian
+! write(*,*) 'dvEdge(200) = ', dvEdge(200)
+
+ do iEdge = 1,nEdges
+ Jacobian = u(k,iEdge)*gradPVn(k,iEdge) + v(k,iEdge)*gradPVt(k,iEdge)
+
+ coef = alpha*m_ke**(-1.5)*m_h*m_Jacobian*dvEdge(iEdge)**3
+ !coef = .5*dt
+
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - coef*Jacobian
+
+ end do
+
+! write(*,*) "alpha = ", alpha
+! write(*,*) "coef = ", coef
+ end do
+ end if
+
+
end subroutine compute_solve_diagnostics
</font>
</pre>