<p><b>xylar@lanl.gov</b> 2010-02-05 16:51:55 -0700 (Fri, 05 Feb 2010)</p><p>* Added global_diagnostics module that computes a few global reductions (more to come!)<br>
* Added Todd's vector reconstruction module that uses radial basis functions<br>
Note: this module requires LAPACK for matrix inversion<br>
* Made changes from the shallow water branch to handle vertexDegree != 3 <br>
* Fixed a bug in fortprintf.c that occasionally put an & on its own line<br>
<br>
-Xylar<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean/Makefile
===================================================================
--- branches/ocean/Makefile 2010-02-02 23:59:45 UTC (rev 113)
+++ branches/ocean/Makefile 2010-02-05 23:51:55 UTC (rev 114)
@@ -33,7 +33,7 @@
CPPINCLUDES = -I../inc -I$(NETCDF)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include
-LIBS = -L$(NETCDF)/lib -lnetcdf
+LIBS = -L$(NETCDF)/lib -lnetcdf -llapack
RM = rm -f
CPP = cpp -C -P -traditional
Modified: branches/ocean/Registry/Registry
===================================================================
--- branches/ocean/Registry/Registry 2010-02-02 23:59:45 UTC (rev 113)
+++ branches/ocean/Registry/Registry 2010-02-05 23:51:55 UTC (rev 114)
@@ -7,6 +7,9 @@
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
namelist real sw_model config_visc 0.0
+namelist character io config_input_name grid.nc
+namelist character io config_output_name output.nc
+namelist character io config_restart_name restart.nc
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
@@ -20,12 +23,13 @@
dim maxEdges2 maxEdges2
dim nVertices nVertices
dim TWO 2
-dim THREE 3
+dim R3 3
+dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
dim nTracers nTracers
#
-# var type name_in_file ( dims ) iro- name_in_code
+# var type name_in_file ( dims ) iro- name_in_code super-array array_class
#
var real xtime ( Time ) ro xtime - -
@@ -66,13 +70,19 @@
var integer cellsOnCell ( maxEdges nCells ) iro cellsOnCell - -
var integer verticesOnCell ( maxEdges nCells ) iro verticesOnCell - -
var integer verticesOnEdge ( TWO nEdges ) iro verticesOnEdge - -
-var integer edgesOnVertex ( THREE nVertices ) iro edgesOnVertex - -
-var integer cellsOnVertex ( THREE nVertices ) iro cellsOnVertex - -
-var real kiteAreasOnVertex ( THREE nVertices ) iro kiteAreasOnVertex - -
+var integer edgesOnVertex ( vertexDegree nVertices ) iro edgesOnVertex - -
+var integer cellsOnVertex ( vertexDegree nVertices ) iro cellsOnVertex - -
+var real kiteAreasOnVertex ( vertexDegree nVertices ) iro kiteAreasOnVertex - -
var real fEdge ( nEdges ) iro fEdge - -
var real fVertex ( nVertices ) iro fVertex - -
var real h_s ( nCells ) iro h_s - -
+# Arrays required for reconstruction of velocity field
+var real matrix_reconstruct ( maxEdges maxEdges nCells ) - matrix_reconstruct - -
+var real normal ( R3 maxEdges nCells ) - normal - -
+var real rbf_value ( maxEdges nCells ) - rbf_value - -
+var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
+
# Boundary conditions: read from input, saved in restart and written to output
var integer uBC ( nVertLevels nEdges ) iro uBC - -
var real u_src ( nVertLevels nEdges ) iro u_src - -
@@ -91,9 +101,27 @@
var real ke ( nVertLevels nCells Time ) o ke - -
var real pv_vertex ( nVertLevels nVertices Time ) o pv_vertex - -
var real pv_cell ( nVertLevels nCells Time ) o pv_cell - -
+var real uReconstructX ( nVertLevels nCells Time ) o uReconstructX - -
+var real uReconstructY ( nVertLevels nCells Time ) o uReconstructY - -
+var real uReconstructZ ( nVertLevels nCells Time ) o uReconstructZ - -
# Other diagnostic variables: neither read nor written to any files
var real vh ( nVertLevels nEdges Time ) - vh - -
var real circulation ( nVertLevels nVertices Time ) - circulation - -
var real gradPVt ( nVertLevels nEdges Time ) - gradPVt - -
var real gradPVn ( nVertLevels nEdges Time ) - gradPVn - -
+
+# xsad 10-02-05:
+# Globally reduced diagnostic variables: only written to output
+var real cellAreaGlobal ( Time ) o cellAreaGlobal - -
+var real edgeAreaGlobal ( Time ) o edgeAreaGlobal - -
+var real triangleAreaGlobal ( Time ) o triangleAreaGlobal - -
+
+var real cellVolumeGlobal ( Time ) o cellVolumeGlobal - -
+var real keGlobal ( Time ) o keGlobal - -
+var real averageLayerThicknessGlobal ( Time ) o averageLayerThicknessGlobal - -
+var real maximumVelocityGlobal ( Time ) o maximumVelocityGlobal - -
+var real minimumLayerThicknessGlobal ( Time ) o minimumLayerThicknessGlobal - -
+var real CFLNumberGlobal ( Time ) o CFLNumberGlobal - -
+# xsad 10-02-05 end
+
Modified: branches/ocean/Registry/fortprintf.c
===================================================================
--- branches/ocean/Registry/fortprintf.c 2010-02-02 23:59:45 UTC (rev 113)
+++ branches/ocean/Registry/fortprintf.c 2010-02-05 23:51:55 UTC (rev 114)
@@ -24,7 +24,11 @@
for(i=0; i<MAX_LINE_LEN-1 && i<nbuf; i++) {
if (fbuffer[i] == '\'' && (fbuffer[i+1] != '\'' || i == nbuf-1)) inquotes = (inquotes + 1) % 2;
if (fbuffer[i] == '</font>
<font color="gray">') nl = i;
- if (fbuffer[i] == ' ') sp = i;
+ // xsad 10-02-03:
+ //if (fbuffer[i] == ' ') sp = i;
+ if ((fbuffer[i] == ' ') && (i != nbuf-1) && (fbuffer[i+1] != '&'))
+ sp = i;
+ // end xsad 10-02-03
}
if (nbuf <= MAX_LINE_LEN) sp = -1;
Modified: branches/ocean/src/Makefile
===================================================================
--- branches/ocean/src/Makefile 2010-02-02 23:59:45 UTC (rev 113)
+++ branches/ocean/src/Makefile 2010-02-05 23:51:55 UTC (rev 114)
@@ -13,6 +13,8 @@
module_dmpar.o \
module_io_input.o \
module_io_output.o \
+ module_global_diagnostics.o \
+ module_vector_reconstruction.o \
module_time_integration.o \
streams.o
@@ -34,10 +36,14 @@
module_sw_solver.o: module_configure.o module_grid_types.o module_io_output.o module_time_integration.o module_dmpar.o module_timer.o
-module_time_integration.o: module_grid_types.o module_configure.o module_dmpar.o
+module_time_integration.o: module_grid_types.o module_configure.o module_dmpar.o module_global_diagnostics.o module_vector_reconstruction.o
module_test_cases.o: module_grid_types.o module_configure.o module_constants.o module_dmpar.o
+module_vector_reconstruction.o: module_grid_types.o module_configure.o module_constants.o
+
+module_global_diagnostics.o: module_grid_types.o module_configure.o module_dmpar.o
+
swmodel: $(OBJS)
$(FC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
Modified: branches/ocean/src/module_dmpar.F
===================================================================
--- branches/ocean/src/module_dmpar.F 2010-02-02 23:59:45 UTC (rev 113)
+++ branches/ocean/src/module_dmpar.F 2010-02-05 23:51:55 UTC (rev 114)
@@ -241,7 +241,102 @@
end subroutine dmpar_sum_int
+ !xsad 10-01-22:
+ subroutine dmpar_sum_real(dminfo, r, rsum)
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ real(kind=RKIND), intent(in) :: r
+ real(kind=RKIND), intent(out) :: rsum
+
+ integer :: mpi_ierr
+
+#ifdef _MPI
+ call MPI_Allreduce(r, rsum, 1, MPI_REALKIND, MPI_SUM, dminfo % comm, mpi_ierr)
+#else
+ rsum = r
+#endif
+
+ end subroutine dmpar_sum_real
+ !xsad 10-01-22 end
+
+ !xsad 10-02-01:
+ subroutine dmpar_min_int(dminfo, i, imin)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: i
+ integer, intent(out) :: imin
+
+ integer :: mpi_ierr
+
+#ifdef _MPI
+ call MPI_Allreduce(i, imin, 1, MPI_INTEGER, MPI_MIN, dminfo % comm, mpi_ierr)
+#else
+ imin = i
+#endif
+
+ end subroutine dmpar_min_int
+
+ subroutine dmpar_min_real(dminfo, r, rmin)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ real(kind=RKIND), intent(in) :: r
+ real(kind=RKIND), intent(out) :: rmin
+
+ integer :: mpi_ierr
+
+#ifdef _MPI
+ call MPI_Allreduce(r, rmin, 1, MPI_REALKIND, MPI_MIN, dminfo % comm, mpi_ierr)
+#else
+ rmin = r
+#endif
+
+ end subroutine dmpar_min_real
+
+ subroutine dmpar_max_int(dminfo, i, imax)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: i
+ integer, intent(out) :: imax
+
+ integer :: mpi_ierr
+
+#ifdef _MPI
+ call MPI_Allreduce(i, imax, 1, MPI_INTEGER, MPI_MAX, dminfo % comm, mpi_ierr)
+#else
+ imax = i
+#endif
+
+ end subroutine dmpar_max_int
+
+ subroutine dmpar_max_real(dminfo, r, rmax)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ real(kind=RKIND), intent(in) :: r
+ real(kind=RKIND), intent(out) :: rmax
+
+ integer :: mpi_ierr
+
+#ifdef _MPI
+ call MPI_Allreduce(r, rmax, 1, MPI_REALKIND, MPI_MAX, dminfo % comm, mpi_ierr)
+#else
+ rmax = r
+#endif
+
+ end subroutine dmpar_max_real
+
+ !xsad 10-02-01 end
+
+
subroutine dmpar_scatter_ints(dminfo, nprocs, noutlist, displs, counts, inlist, outlist)
implicit none
Added: branches/ocean/src/module_global_diagnostics.F
===================================================================
--- branches/ocean/src/module_global_diagnostics.F (rev 0)
+++ branches/ocean/src/module_global_diagnostics.F 2010-02-05 23:51:55 UTC (rev 114)
@@ -0,0 +1,180 @@
+module global_diagnostics
+
+ use grid_types
+ use configure
+ use constants
+ use dmpar
+
+ implicit none
+ save
+ public
+
+ contains
+
+ subroutine computeGlobalDiagnostics(dminfo, state, grid, dt)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ type (grid_state), intent(inout) :: state
+ type (grid_meta), intent(in) :: grid
+ real (kind=RKIND), intent(in) :: dt
+
+ type (block_type), pointer :: block
+ integer :: nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve
+ real (kind=RKIND) :: cellAreaGlobal, edgeAreaGlobal, triangleAreaGlobal
+ real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle
+ real (kind=RKIND), dimension(:,:), pointer :: h, ke, u
+ real (kind=RKIND) :: cellVolumeGlobal, keGlobal, averageLayerThicknessGlobal, maximumVelocityGlobal, &
+ minimumLayerThicknessGlobal, CFLNumberGlobal
+ real (kind=RKIND) :: localCFL, localSum
+ integer :: elementIndex
+ integer :: timeLevel
+
+ nVertLevels = grid % nVertLevels
+ nCellsSolve = grid % nCellsSolve
+ nEdgesSolve = grid % nEdgesSolve
+ nVerticesSolve = grid % nVerticesSolve
+
+ areaCell => grid % areaCell % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaTriangle => grid % areaTriangle % array
+
+ h => state % h % array
+ ke => state % ke % array
+ u => state % u % array
+
+ localSum = sum(areaCell(1:nCellsSolve))
+ call dmpar_sum_real(dminfo, localSum, cellAreaGlobal)
+
+ localSum = sum(dcEdge(1:nEdgesSolve)*dvEdge(1:nEdgesSolve))
+ call dmpar_sum_real(dminfo, localSum, edgeAreaGlobal)
+
+ localSum = sum(areaTriangle(1:nVerticesSolve))
+ call dmpar_sum_real(dminfo, localSum, triangleAreaGlobal)
+
+ call computeAreaWeightedGlobalSum(dminfo, nVertLevels, nCellsSolve, &
+ areaCell(1:nCellsSolve), h(:,1:nCellsSolve), cellVolumeGlobal)
+
+ averageLayerThicknessGlobal = cellVolumeGlobal/(cellAreaGlobal*nVertLevels)
+
+ call computeVolumeWeightedGlobalSum(dminfo, nVertLevels, nCellsSolve, &
+ areaCell(1:nCellsSolve), h(:,1:nCellsSolve), ke(:,1:nCellsSolve), keGlobal)
+
+ call computeGlobalMax(dminfo, nVertLevels, nEdgesSolve, &
+ u(:,1:nEdgesSolve), maximumVelocityGlobal)
+
+ call computeGlobalMin(dminfo, nVertLevels, nCellsSolve, &
+ h(:,1:nEdgesSolve), minimumLayerThicknessGlobal)
+
+ localCFL = 0.0
+ do elementIndex = 1,nEdgesSolve
+ localCFL = max(localCFL, maxval(dt*u(:,elementIndex)/dcEdge(elementIndex)))
+ end do
+ call dmpar_max_real(dminfo, localCFL, CFLNumberGlobal)
+
+
+ state % cellAreaGlobal % scalar = cellVolumeGlobal
+ state % edgeAreaGlobal % scalar = edgeAreaGlobal
+ state % triangleAreaGlobal % scalar = triangleAreaGlobal
+
+ state % cellVolumeGlobal % scalar = cellVolumeGlobal
+ state % keGlobal % scalar = keGlobal
+ state % averageLayerThicknessGlobal % scalar = averageLayerThicknessGlobal
+ state % maximumVelocityGlobal % scalar = maximumVelocityGlobal
+ state % minimumLayerThicknessGlobal % scalar = minimumLayerThicknessGlobal
+ state % CFLNumberGlobal % scalar = CFLNumberGlobal
+
+ end subroutine computeGlobalDiagnostics
+
+ subroutine computeGlobalSum(dminfo, nVertLevels, nElements, field, globalSum)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalSum
+
+ real (kind=RKIND) :: localSum
+
+ localSum = sum(field)
+ call dmpar_sum_real(dminfo, localSum, globalSum)
+
+ end subroutine computeGlobalSum
+
+ subroutine computeAreaWeightedGlobalSum(dminfo, nVertLevels, nElements, areas, field, globalSum)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nElements), intent(in) :: areas
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalSum
+
+ integer :: elementIndex
+ real (kind=RKIND) :: localSum
+
+ localSum = 0.
+ do elementIndex = 1, nElements
+ localSum = localSum + areas(elementIndex) * sum(field(:,elementIndex))
+ end do
+
+ call dmpar_sum_real(dminfo, localSum, globalSum)
+
+ end subroutine computeAreaWeightedGlobalSum
+
+ subroutine computeVolumeWeightedGlobalSum(dminfo, nVertLevels, nElements, areas, h, field, globalSum)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nElements), intent(in) :: areas
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: h
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalSum
+
+ real (kind=RKIND), dimension(nVertLevels, nElements) :: hTimesField
+
+ hTimesField = h*field
+
+ call computeAreaWeightedGlobalSum(dminfo, nVertLevels, nElements, areas, hTimesField, globalSum)
+
+ end subroutine computeVolumeWeightedGlobalSum
+
+ subroutine computeGlobalMin(dminfo, nVertLevels, nElements, field, globalMin)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMin
+
+ real (kind=RKIND) :: localMin
+
+ localMin = minval(field)
+ call dmpar_min_real(dminfo, localMin, globalMin)
+
+ end subroutine computeGlobalMin
+
+ subroutine computeGlobalMax(dminfo, nVertLevels, nElements, field, globalMax)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMax
+
+ real (kind=RKIND) :: localMax
+
+ localMax = maxval(field)
+ call dmpar_max_real(dminfo, localMax, globalMax)
+
+ end subroutine computeGlobalMax
+
+end module global_diagnostics
Modified: branches/ocean/src/module_io_input.F
===================================================================
--- branches/ocean/src/module_io_input.F 2010-02-02 23:59:45 UTC (rev 113)
+++ branches/ocean/src/module_io_input.F 2010-02-05 23:51:55 UTC (rev 114)
@@ -91,9 +91,9 @@
integer :: idxTdiff
if (config_do_restart) then
- input_obj % filename = 'restart.nc'
+ input_obj % filename = trim(config_restart_name)
else
- input_obj % filename = 'grid.nc'
+ input_obj % filename = trim(config_input_name)
end if
call io_input_init(input_obj, domain % dminfo)
@@ -207,9 +207,9 @@
cellsOnVertexField % ioinfo % fieldName = 'cellsOnVertex'
cellsOnVertexField % ioinfo % start(1) = 1
cellsOnVertexField % ioinfo % start(2) = readVertexStart
- cellsOnVertexField % ioinfo % count(1) = 3
+ cellsOnVertexField % ioinfo % count(1) = vertexDegree
cellsOnVertexField % ioinfo % count(2) = nReadVertices
- allocate(cellsOnVertexField % array(3,nReadVertices))
+ allocate(cellsOnVertexField % array(vertexDegree,nReadVertices))
call io_input_field(input_obj, cellsOnVertexField)
@@ -397,14 +397,14 @@
! Work out which edges and vertices are owned by this process, and which are ghost
!
allocate(cellsOnEdge_2Halo(2,nlocal_edges))
- allocate(cellsOnVertex_2Halo(3,nlocal_vertices))
+ allocate(cellsOnVertex_2Halo(vertexDegree,nlocal_vertices))
call dmpar_alltoall_field(domain % dminfo, cellsOnEdgeField % array, cellsOnEdge_2Halo, &
2, size(cellsOnEdgeField % array, 2), nlocal_edges, &
sendEdgeList, recvEdgeList)
call dmpar_alltoall_field(domain % dminfo, cellsOnVertexField % array, cellsOnVertex_2Halo, &
- 3, size(cellsOnVertexField % array, 2), nlocal_vertices, &
+ vertexDegree, size(cellsOnVertexField % array, 2), nlocal_vertices, &
sendVertexList, recvVertexList)
@@ -413,7 +413,7 @@
2, nlocal_edges, cellsOnEdge_2Halo, local_edge_list, ghostEdgeStart)
call block_decomp_partitioned_edge_list(block_graph_2Halo % nVertices, &
block_graph_2Halo % vertexID(1:block_graph_2Halo % nVertices), &
- 3, nlocal_vertices, cellsOnVertex_2Halo, local_vertex_list, ghostVertexStart)
+ vertexDegree, nlocal_vertices, cellsOnVertex_2Halo, local_vertex_list, ghostVertexStart)
! At this point, local_edge_list(1;ghostEdgeStart-1) contains all of the owned edges for this block
@@ -501,7 +501,7 @@
!
! If doing a restart, we need to decide which time slice to read from the
- ! restart.nc file
+ ! restart file
!
if (input_obj % rdLocalTime <= 0) then
write(0,*) 'Error: Couldn''t find any times in restart file.'
@@ -649,7 +649,7 @@
end do
do i=1,domain % blocklist % mesh % nVertices
- do j=1,3
+ do j=1,vertexDegree
k = binary_search(cellIDSorted, 2, 1, domain % blocklist % mesh % nCells, &
domain % blocklist % mesh % cellsOnVertex % array(j,i))
Modified: branches/ocean/src/module_io_output.F
===================================================================
--- branches/ocean/src/module_io_output.F 2010-02-02 23:59:45 UTC (rev 113)
+++ branches/ocean/src/module_io_output.F 2010-02-05 23:51:55 UTC (rev 114)
@@ -3,6 +3,7 @@
use grid_types
use dmpar
use sort
+ use configure
integer, parameter :: OUTPUT = 1
integer, parameter :: RESTART = 2
@@ -74,10 +75,10 @@
nVertLevelsGlobal = block_ptr % mesh % nVertLevels
if (trim(stream) == 'OUTPUT') then
- output_obj % filename = 'output.nc'
+ output_obj % filename = trim(config_output_name)
output_obj % stream = OUTPUT
else if (trim(stream) == 'RESTART') then
- output_obj % filename = 'restart.nc'
+ output_obj % filename = trim(config_restart_name)
output_obj % stream = RESTART
end if
@@ -149,8 +150,8 @@
allocate(cellsOnEdge(2, domain % blocklist % mesh % nEdgesSolve))
allocate(verticesOnEdge(2, domain % blocklist % mesh % nEdgesSolve))
allocate(edgesOnEdge(2 * domain % blocklist % mesh % maxEdges, domain % blocklist % mesh % nEdgesSolve))
- allocate(cellsOnVertex(3, domain % blocklist % mesh % nVerticesSolve))
- allocate(edgesOnVertex(3, domain % blocklist % mesh % nVerticesSolve))
+ allocate(cellsOnVertex(domain % blocklist % mesh % vertexDegree, domain % blocklist % mesh % nVerticesSolve))
+ allocate(edgesOnVertex(domain % blocklist % mesh % vertexDegree, domain % blocklist % mesh % nVerticesSolve))
!
@@ -191,7 +192,7 @@
end do
end do
do i=1,domain % blocklist % mesh % nVerticesSolve
- do j=1,3
+ do j=1,domain % blocklist % mesh % vertexDegree
cellsOnVertex(j,i) = domain % blocklist % mesh % indexToCellID % array( &
domain % blocklist % mesh % cellsOnVertex % array(j,i))
edgesOnVertex(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
Modified: branches/ocean/src/module_time_integration.F
===================================================================
--- branches/ocean/src/module_time_integration.F 2010-02-02 23:59:45 UTC (rev 113)
+++ branches/ocean/src/module_time_integration.F 2010-02-05 23:51:55 UTC (rev 114)
@@ -4,6 +4,10 @@
use configure
use constants
use dmpar
+ ! xsad 10-02-05:
+ use global_diagnostics
+ use vector_reconstruction
+ ! xsad 10-02-05 end
contains
@@ -213,6 +217,12 @@
call compute_solve_diagnostics(dt, block % time_levs(2) % state, block % mesh)
+ call reconstruct(block % time_levs(2) % state, block % mesh)
+
+ ! xsad 10-02-05:
+ call computeGlobalDiagnostics(domain % dminfo, block % time_levs(2) % state, block % mesh, dt)
+ ! xsad 10-02-05 end
+
block => block % next
end do
@@ -620,22 +630,20 @@
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
- do iVertex = 1,nVertices
- if (cellsOnVertex(1,iVertex) > 0 .and. &
- cellsOnVertex(2,iVertex) > 0 .and. &
- cellsOnVertex(3,iVertex) > 0 &
- ) then
- do k=1,nVertLevels
- h_vertex = 0.0
- do i=1,3
- h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
- end do
- h_vertex = h_vertex / areaTriangle(iVertex)
-
- pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
+ VTX_LOOP: do iVertex = 1,nVertices
+ do i=1,grid % vertexDegree
+ if (cellsOnVertex(i,iVertex) <= 0) cycle VTX_LOOP
+ end do
+ do k=1,nVertLevels
+ h_vertex = 0.0
+ do i=1,grid % vertexDegree
+ h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
- end if
- end do
+ h_vertex = h_vertex / areaTriangle(iVertex)
+
+ pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
+ end do
+ end do VTX_LOOP
! tdr
@@ -658,7 +666,7 @@
!
pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
- do i=1,3
+ do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
if(iEdge > 0) then
do k=1,nVertLevels
@@ -687,7 +695,7 @@
!
pv_cell(:,:) = 0.0
do iVertex = 1, nVertices
- do i=1,3
+ do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
if( iCell > 0) then
do k = 1,nVertLevels
@@ -786,5 +794,4 @@
end subroutine enforce_uBC
-
end module time_integration
Added: branches/ocean/src/module_vector_reconstruction.F
===================================================================
--- branches/ocean/src/module_vector_reconstruction.F (rev 0)
+++ branches/ocean/src/module_vector_reconstruction.F 2010-02-05 23:51:55 UTC (rev 114)
@@ -0,0 +1,384 @@
+ module vector_reconstruction
+
+ use grid_types
+ use configure
+ use constants
+
+ implicit none
+
+ public :: init_reconstruct, reconstruct
+
+ contains
+
+ subroutine init_reconstruct(grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: reconstruct vector field at cell centers based on radial basis functions
+ !
+ ! Input: grid meta data and vector component data residing at cell edges
+ !
+ ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_meta), intent(inout) :: grid
+
+ ! temporary arrays needed in the (to be constructed) init procedure
+ integer :: nCells, nEdges, nVertLevels, nCellsSolve
+ integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, cellsOnEdge
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer :: iEdge, iCell, k, cell1, cell2, EdgeMax, j, i, npoints
+ integer :: lwork, info
+ integer, allocatable, dimension(:) :: ipvt
+ real (kind=RKIND), allocatable, dimension(:) :: work
+ real (kind=RKIND), dimension(:), pointer :: dcEdge, xCell, yCell, zCell
+ real (kind=RKIND) :: r, t, v, X1(3), X2(3), alpha
+ real (kind=RKIND), allocatable, dimension(:,:,:) :: xLoc
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: matrix_reconstruct
+ real (kind=RKIND), dimension(:,:,:), pointer :: normal
+ real (kind=RKIND), dimension(:,:), pointer :: rbf_value
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+ real (kind=RKIND), allocatable, dimension(:,:) :: mwork
+
+ !========================================================
+ ! arrays filled and saved during init procedure
+ !========================================================
+ matrix_reconstruct => grid % matrix_reconstruct % array
+ normal => grid % normal % array
+ rbf_value => grid % rbf_value % array
+ coeffs_reconstruct => grid % coeffs_reconstruct % array
+
+ !========================================================
+ ! temporary variables needed for init procedure
+ !========================================================
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ cellsOnCell => grid % cellsOnCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ dcEdge => grid % dcEdge % array
+ nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+
+ ! allocate arrays
+ EdgeMax = maxval(nEdgesOnCell)
+ allocate(xLoc(3,EdgeMax,nCells))
+ allocate(work(EdgeMax*EdgeMax))
+ allocate(ipvt(EdgeMax))
+
+ ! init arrays
+ matrix_reconstruct = 0.0
+ normal = 0.0
+ rbf_value = 0.0
+
+ ! loop over all cells to be solved on this block
+ do iCell=1,nCellsSolve
+
+ npoints = nEdgesOnCell(iCell) ! only loop over the number of edges for this cell
+
+ ! fill normal vector and xLoc arrays
+ ! X1 is the location of the reconstruction, X2 are the neighboring centers, xLoc is the edge positions
+ cell1 = iCell
+ X1(1) = xCell(cell1)
+ X1(2) = yCell(cell1)
+ X1(3) = zCell(cell1)
+ do j=1,npoints
+ cell2 = cellsOnCell(j,iCell)
+ iEdge = edgesOnCell(j,iCell)
+ X2(1) = xCell(cell2)
+ X2(2) = yCell(cell2)
+ X2(3) = zCell(cell2)
+ if(cell2.gt.cell1) then
+ normal(:,j,iCell) = X2(:) - X1(:)
+ else
+ normal(:,j,iCell) = X1(:) - X2(:)
+ endif
+ call unit_vector_in_R3(normal(:,j,iCell))
+ xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
+ enddo
+
+ alpha = 0.0
+ do i=1,npoints
+ call get_distance(xLoc(:,i,iCell),X1(:),r)
+ alpha = alpha + r
+ enddo
+ alpha = 4*alpha / npoints
+ do j=1,npoints
+ do i=1,npoints
+ call get_distance(xLoc(:,i,iCell),xLoc(:,j,iCell),r)
+ r = r / alpha
+ call evaluate_rbf(r,t)
+ call get_dotproduct(normal(:,i,iCell),normal(:,j,iCell),v)
+ matrix_reconstruct(i,j,iCell) = v*t
+ enddo
+ enddo
+
+ ! save value of RBF when evaluated between reconstruction location and edge locations
+ do j=1,npoints
+ call get_distance(xLoc(:,j,iCell), X1(:), r)
+ r = r / alpha
+ call evaluate_rbf(r,t)
+ rbf_value(j,iCell) = t
+ enddo
+
+ ! invert matrix_reconstruct using lapack routines
+ allocate(mwork(npoints,npoints))
+ lwork = npoints*npoints
+ mwork(:,:) = matrix_reconstruct(1:npoints,1:npoints,iCell)
+ call dgetrf(npoints, npoints, mwork, npoints, ipvt, info)
+ if(info.ne.0) then
+ write(6,*) info, 'dgetrf'
+ stop
+ endif
+ call dgetri(npoints, mwork, npoints, ipvt, work, lwork, info)
+ matrix_reconstruct(1:npoints,1:npoints,iCell) = mwork(:,:)
+ if(info.ne.0) then
+ write(6,*) info, 'dgetri'
+ stop
+ endif
+ deallocate(mwork)
+
+ ! compute the coefficients for reconstructing uX, uY, uZ at cell centers from u_i normal to edges
+ ! uX = sum_j(coeffs(1,j) * u_j) (similarly for Y,Z)
+ ! coeffs(:,j) = sum_i(rbf_values(i) * normal(:,i) * matrix(i,j))
+ do j=1,npoints
+ coeffs_reconstruct(:,j,iCell) = 0.0
+ do i=1,npoints
+ coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) + rbf_value(i,iCell) * normal(:,i,iCell) &
+ * matrix_reconstruct(i,j,iCell)
+ enddo
+ enddo
+
+ enddo ! iCell
+
+ deallocate(ipvt)
+ deallocate(work)
+ deallocate(xLoc)
+
+ end subroutine init_reconstruct
+
+
+ subroutine reconstruct(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: reconstruct vector field at cell centers based on radial basis functions
+ !
+ ! Input: grid meta data and vector component data residing at cell edges
+ !
+ ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_state), intent(inout) :: s
+ type (grid_meta), intent(in) :: grid
+
+ ! temporary arrays needed in the compute procedure
+ integer :: nCells, nEdges, nVertLevels, nCellsSolve
+ integer, dimension(:,:), pointer :: edgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+ integer :: iCell,iEdge, i, j, k, npoints, EdgeMax
+ real (kind=RKIND) :: r, t, p!, t1(3), t2(3), t3(3), east(3)
+ real (kind=RKIND), dimension(:,:), pointer :: u
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
+ real (kind=RKIND), allocatable, dimension(:,:) :: rhs, c
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: matrix_reconstruct
+ real (kind=RKIND), dimension(:,:,:), pointer :: normal
+ real (kind=RKIND), dimension(:,:), pointer :: rbf_value
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+
+ ! stored arrays used during compute procedure
+ matrix_reconstruct => grid % matrix_reconstruct % array
+ normal => grid % normal % array
+ rbf_value => grid % rbf_value % array
+ coeffs_reconstruct => grid % coeffs_reconstruct % array
+
+ ! temporary variables
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+
+ u => s % u % array
+ uReconstructX => s % uReconstructX % array
+ uReconstructY => s % uReconstructY % array
+ uReconstructZ => s % uReconstructZ % array
+
+ ! allocate space for temporary arrays
+ EdgeMax = maxval(nEdgesOnCell)
+ allocate(c(nVertLevels,EdgeMax))
+ allocate(rhs(nVertLevels,EdgeMax))
+
+ ! init the intent(out)
+ uReconstructX = 0.0
+ uReconstructY = 0.0
+ uReconstructZ = 0.0
+
+ ! loop over cell centers
+ do iCell=1,nCellsSolve
+
+! t2(1) = xCell(iCell)
+! t2(2) = yCell(iCell)
+! t2(3) = zCell(iCell)
+! t3(1) = 0
+! t3(2) = 0
+! t3(3) = 1.0
+! call cross_product_in_R3(t3(:), t2(:), east(:))
+! call unit_vector_in_R3(east(:))
+
+ npoints = nEdgesOnCell(iCell)
+
+ ! fill the RHS of the matrix equation
+ ! for testing purposes, fill rhs with test function
+ ! test function assumes uniform flow to the east
+ rhs = 0.0
+ do j=1,npoints
+ iEdge = edgesOnCell(j,iCell)
+ do k=1,nVertLevels
+ rhs(k,j) = u(k,iEdge)
+ ! call get_dotproduct(normal(:,j,iCell),east(:),rhs(k,j))
+ enddo
+ enddo
+
+ ! compute c by multiplying matrix_reconstruct * rhs (Eq 8 in Tex document)
+ c = 0.0
+ do i=1,npoints
+ do j=1,npoints
+ do k=1,nVertLevels
+ c(k,i) = c(k,i) + matrix_reconstruct(i,j,iCell)*rhs(k,j)
+ enddo
+ enddo
+ enddo
+
+ ! reconstruct the velocity field at point X1 (Eq 6 in Tex document)
+! do i=1,npoints
+! do k=1,nVertLevels
+! uReconstructX(k,iCell) = uReconstructX(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(1,i,iCell)
+! uReconstructY(k,iCell) = uReconstructY(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(2,i,iCell)
+! uReconstructZ(k,iCell) = uReconstructZ(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(3,i,iCell)
+! enddo
+! enddo
+!
+! t1(1) = uReconstructX(1,iCell)
+! t1(2) = uReconstructY(1,iCell)
+! t1(3) = uReconstructZ(1,iCell)
+! call get_dotproduct(east, t1, p)
+! write(6,10) iCell, p
+! write(6,10) iCell, t1(:)
+! write(6,10) iCell, east(:)
+! write(6,*)
+! 10 format(i6, 3e15.5)
+
+ ! reconstruct the velocity field at point X1 (Eq 6 in Tex document)
+! uReconstructX(:,iCell) = 0.0
+! uReconstructY(:,iCell) = 0.0
+! uReconstructZ(:,iCell) = 0.0
+ do i=1,npoints
+ uReconstructX(:,iCell) = uReconstructX(:,iCell) + coeffs_reconstruct(1,i,iCell) * rhs(:,i)
+ uReconstructY(:,iCell) = uReconstructY(:,iCell) + coeffs_reconstruct(2,i,iCell) * rhs(:,i)
+ uReconstructZ(:,iCell) = uReconstructZ(:,iCell) + coeffs_reconstruct(3,i,iCell) * rhs(:,i)
+ enddo
+
+! t1(1) = uReconstructX(1,iCell)
+! t1(2) = uReconstructY(1,iCell)
+! t1(3) = uReconstructZ(1,iCell)
+! call get_dotproduct(east, t1, p)
+! write(6,10) iCell, p
+! write(6,10) iCell, t1(:)
+! write(6,10) iCell, east(:)
+! write(6,*)
+
+
+ enddo ! iCell
+
+ ! deallocate
+ deallocate(rhs)
+ deallocate(c)
+
+! stop
+
+ end subroutine reconstruct
+
+ subroutine get_distance(x1,x2,xout)
+ implicit none
+ real(kind=RKIND), intent(in) :: x1(3), x2(3)
+ real(kind=RKIND), intent(out) :: xout
+ xout = sqrt( (x1(1)-x2(1))**2 + (x1(2)-x2(2))**2 + (x1(3)-x2(3))**2 )
+ end subroutine get_distance
+
+ subroutine get_dotproduct(x1,x2,xout)
+ implicit none
+ real(kind=RKIND), intent(in) :: x1(3), x2(3)
+ real(kind=RKIND), intent(out) :: xout
+ xout = x1(1)*x2(1) + x1(2)*x2(2) + x1(3)*x2(3)
+ end subroutine get_dotproduct
+
+
+ subroutine unit_vector_in_R3(xin)
+ implicit none
+ real (kind=RKIND), intent(inout) :: xin(3)
+ real (kind=RKIND) :: mag
+ mag = sqrt(xin(1)**2+xin(2)**2+xin(3)**2)
+ xin(:) = xin(:) / mag
+ end subroutine unit_vector_in_R3
+
+
+ subroutine evaluate_rbf(r,xout)
+ real(kind=RKIND), intent(in) :: r
+ real(kind=RKIND), intent(out) :: xout
+
+ ! Gaussian
+ xout = exp(-r**2)
+
+ ! multiquadrics
+ ! xout = 1.0 / sqrt(1.0**2 + r**2)
+
+ ! other
+ ! xout = 1.0 / (1.0**2 + r**2)
+
+ end subroutine evaluate_rbf
+
+!======================================================================
+! BEGINNING OF CROSS_PRODUCT_IN_R3
+!======================================================================
+ subroutine cross_product_in_R3(p_1,p_2,p_out)
+
+!-----------------------------------------------------------------------
+! PURPOSE: compute p_1 cross p_2 and place in p_out
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(in)
+!-----------------------------------------------------------------------
+ real (kind=RKIND), intent(in) :: &
+ p_1 (3), &
+ p_2 (3)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+ real (kind=RKIND), intent(out) :: &
+ p_out (3)
+
+ p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
+ p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
+ p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
+
+ end subroutine cross_product_in_R3
+!======================================================================
+! END OF CROSS_PRODUCT_IN_R3
+!======================================================================
+
+
+end module vector_reconstruction
</font>
</pre>