<p><b>ringler@lanl.gov</b> 2010-01-28 21:32:22 -0700 (Thu, 28 Jan 2010)</p><p><br>
adding a vector reconstruction module to SWM<br>
<br>
the reconstruction is based on radial basis functions<br>
<br>
somewhat surprisingly, the module appears to be returning approximately correct results when compared to a test function.<br>
<br>
please note: currently the module is using a test function on the RHS<br>
the module needs to be separted into init and compute procedures<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/Makefile
===================================================================
--- trunk/swmodel/Makefile 2010-01-07 21:57:54 UTC (rev 111)
+++ trunk/swmodel/Makefile 2010-01-29 04:32:22 UTC (rev 112)
@@ -50,7 +50,7 @@
"SCC = gcc" \
"FFLAGS = -real-size 64 -O3" \
"CFLAGS = -O3 -m64" \
- "LDFLAGS = -O3" \
+ "LDFLAGS = -O3 -L$MKLPATH -I$MKLINCLUDE -I$MKLPATH -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
gfortran:
Modified: trunk/swmodel/Registry/Registry
===================================================================
--- trunk/swmodel/Registry/Registry 2010-01-07 21:57:54 UTC (rev 111)
+++ trunk/swmodel/Registry/Registry 2010-01-29 04:32:22 UTC (rev 112)
@@ -93,6 +93,9 @@
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 - -
Modified: trunk/swmodel/src/Makefile
===================================================================
--- trunk/swmodel/src/Makefile 2010-01-07 21:57:54 UTC (rev 111)
+++ trunk/swmodel/src/Makefile 2010-01-29 04:32:22 UTC (rev 112)
@@ -19,12 +19,13 @@
module_io_input.o \
module_io_output.o \
module_time_integration.o \
+ module_vector_reconstruction.o \
$(ZOLTANOBJ) \
streams.o
all: swmodel
-swmodel.o: module_configure.o module_dmpar.o module_grid_types.o module_test_cases.o module_io_input.o module_sw_solver.o module_timer.o
+swmodel.o: module_configure.o module_dmpar.o module_grid_types.o module_test_cases.o module_io_input.o module_sw_solver.o module_timer.o
module_configure.o: module_dmpar.o
@@ -40,10 +41,12 @@
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_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
+
swmodel: $(OBJS)
$(FC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
Modified: trunk/swmodel/src/module_time_integration.F
===================================================================
--- trunk/swmodel/src/module_time_integration.F 2010-01-07 21:57:54 UTC (rev 111)
+++ trunk/swmodel/src/module_time_integration.F 2010-01-29 04:32:22 UTC (rev 112)
@@ -1,5 +1,6 @@
module time_integration
+ use vector_reconstruction
use grid_types
use configure
use constants
@@ -213,6 +214,8 @@
call compute_solve_diagnostics(dt, block % time_levs(2) % state, block % mesh)
+ call reconstruct(block % time_levs(2) % state, block % mesh)
+
block => block % next
end do
Added: trunk/swmodel/src/module_vector_reconstruction.F
===================================================================
--- trunk/swmodel/src/module_vector_reconstruction.F (rev 0)
+++ trunk/swmodel/src/module_vector_reconstruction.F 2010-01-29 04:32:22 UTC (rev 112)
@@ -0,0 +1,276 @@
+ module vector_reconstruction
+
+ use grid_types
+ use configure
+ use constants
+
+ public :: reconstruct
+
+ contains
+
+ 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 ! s will be intent(inout) for the compute procedure
+ type (grid_meta), intent(in) :: grid ! grid will be intent(in) for the init procedure
+
+ ! arrays needed to be initialized and saved in the (to be constructed) init procedure
+ real (kind=RKIND), allocatable, dimension(:,:,:) :: matrix_reconstruct
+ real (kind=RKIND), allocatable, dimension(:,:,:) :: normal
+ real (kind=RKIND), allocatable, dimension(:,:) :: rbf_value, mwork
+
+ ! 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
+
+ ! temporary arrays needed in the compute procedure (some are commented to prevent double declarations)
+! integer :: nCells, nEdges, nVertLevels
+! integer, dimension(:,:), pointer :: edgesOnCell
+! integer, dimension(:), pointer :: nEdgesOnCell
+! integer :: iCell,iEdge, i, j, k, npoints
+ real (kind=RKIND), dimension(:,:), pointer :: u
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
+ real (kind=RKIND), allocatable, dimension(:,:) :: rhs, c
+
+ !========================================================
+ ! 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
+
+
+ !========================================
+ ! needed for compute procedure
+ !========================================
+ u => s % u % array
+ uReconstructX => s % uReconstructX % array
+ uReconstructY => s % uReconstructY % array
+ uReconstructZ => s % uReconstructZ % array
+
+
+ ! allocate arrays, all allocations occur in init procedure
+ EdgeMax = maxval(nEdgesOnCell)
+ allocate(normal(3,EdgeMax,nCells))
+ allocate(matrix_reconstruct(EdgeMax, EdgeMax, nCells))
+ allocate(xLoc(3,EdgeMax,nCells))
+ allocate(rbf_value(EdgeMax, nCells))
+ allocate(work(EdgeMax*EdgeMax))
+ allocate(rhs(nVertLevels,EdgeMax))
+ allocate(c(nVertLevels,EdgeMax))
+ allocate(ipvt(EdgeMax))
+
+
+ ! fill normal vector and xLoc arrays
+ ! X1 is the location of the reconstruction, X2 are the neighboring centers, xLoc is the edge positions
+ do iCell=1,nCellsSolve
+ cell1 = iCell
+ X1(1) = xCell(cell1)
+ X1(2) = yCell(cell1)
+ X1(3) = zCell(cell1)
+ do j=1,nEdgesOnCell(iCell)
+ 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
+ enddo
+
+ ! fill matrix_reconstruct array
+ matrix_reconstruct = 0.0
+ rbf_value = 0.0
+ do iCell=1,nCellsSolve
+ cell1 = iCell
+ X1(1) = xCell(cell1)
+ X1(2) = yCell(cell1)
+ X1(3) = zCell(cell1)
+ npoints = nEdgesOnCell(iCell) ! only loop over the number of edges for this cell
+ alpha = 0.0
+ do i=1,npoints
+ call get_distance(xLoc(:,i,iCell),X1(:),r)
+ alpha = alpha + r
+ enddo
+ alpha = 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)
+
+ enddo ! iCell
+
+!=================================================================
+! this ends the init procedure, the compute procedure begins below
+!=================================================================
+
+ ! init the intent(out)
+ uReconstructX = 0.0
+ uReconstructY = 0.0
+ uReconstructZ = 0.0
+
+ ! loop over cell centers
+ do iCell=1,nCellsSolve
+ npoints = nEdgesOnCell(iCell)
+
+ ! fill the RHS of the matrix equation
+ ! for testing purposes, fill rhs with test function
+ ! test function assumes uniform flow in direction of normal vector of edge 1
+ do j=1,npoints
+ iEdge = edgesOnCell(j,iCell)
+ do k=1,nVertLevels
+ ! rhs(k,j) = u(k,iEdge)
+ r = normal(1,1,iCell)*normal(1,j,iCell)+normal(2,1,iCell)*normal(2,j,iCell)+normal(3,1,iCell)*normal(3,j,iCell)
+ rhs(k,j) = r
+ 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
+
+ ! recontruct 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
+
+ ! uncomment following lines to compare solution to test function
+ ! write(6,10) iCell, uReconstructX(1,iCell), uReconstructY(1,iCell), uReconstructZ(1,iCell)
+ ! write(6,10) iCell, normal(1,1,iCell), normal(2,1,iCell), normal(3,1,iCell)
+ ! write(6,*)
+ !10 format(i6,3f16.6)
+
+ enddo ! iCell
+
+ ! deallocate (put at end of init procedure)
+ deallocate(xLoc)
+ deallocate(work)
+ deallocate(ipvt)
+
+ ! deallocate (put at end of compute procedure)
+ deallocate(rhs)
+ deallocate(c)
+
+ ! deallocate (these deallocates will be removed when init procedure is implemented)
+ deallocate(normal)
+ deallocate(matrix_reconstruct)
+ deallocate(rbf_value)
+
+ contains
+
+ 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(xin,xout)
+ real(kind=RKIND), intent(in) :: xin
+ 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
+
+
+
+ end subroutine reconstruct
+
+
+end module vector_reconstruction
</font>
</pre>