<p><b>ringler@lanl.gov</b> 2010-02-02 16:59:45 -0700 (Tue, 02 Feb 2010)</p><p><br>
separated the vector reconstruction into an init and compute procedure.<br>
<br>
more testing is required before the reconstruction can be used as a part of the solution procedure (i.e. for the physics)<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/Registry/Registry
===================================================================
--- trunk/swmodel/Registry/Registry 2010-01-29 04:32:22 UTC (rev 112)
+++ trunk/swmodel/Registry/Registry 2010-02-02 23:59:45 UTC (rev 113)
@@ -23,6 +23,7 @@
dim maxEdges2 maxEdges2
dim nVertices nVertices
dim TWO 2
+dim R3 3
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
dim nTracers nTracers
@@ -76,6 +77,11 @@
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 - -
+
# Boundary conditions: read from input, saved in restart and written to output
var integer uBC ( nVertLevels nEdges ) iro uBC - -
Modified: trunk/swmodel/src/Makefile
===================================================================
--- trunk/swmodel/src/Makefile 2010-01-29 04:32:22 UTC (rev 112)
+++ trunk/swmodel/src/Makefile 2010-02-02 23:59:45 UTC (rev 113)
@@ -39,7 +39,7 @@
module_io_output.o: module_grid_types.o module_dmpar.o module_sort.o module_configure.o
-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_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_vector_reconstruction.o
module_time_integration.o: module_grid_types.o module_configure.o module_dmpar.o module_vector_reconstruction.o
Modified: trunk/swmodel/src/module_sw_solver.F
===================================================================
--- trunk/swmodel/src/module_sw_solver.F 2010-01-29 04:32:22 UTC (rev 112)
+++ trunk/swmodel/src/module_sw_solver.F 2010-02-02 23:59:45 UTC (rev 113)
@@ -6,6 +6,7 @@
use configure
use dmpar
use timer
+ use vector_reconstruction
integer :: output_frame
integer :: restart_frame
@@ -44,6 +45,7 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
call compute_solve_diagnostics(dt,block_ptr % time_levs(1) % state, block_ptr % mesh)
+ call init_reconstruct(block_ptr % mesh)
if (.not. config_do_restart) block_ptr % time_levs(1) % state % xtime % scalar = 0.0
block_ptr => block_ptr % next
end do
Modified: trunk/swmodel/src/module_vector_reconstruction.F
===================================================================
--- trunk/swmodel/src/module_vector_reconstruction.F 2010-01-29 04:32:22 UTC (rev 112)
+++ trunk/swmodel/src/module_vector_reconstruction.F 2010-02-02 23:59:45 UTC (rev 113)
@@ -4,11 +4,11 @@
use configure
use constants
- public :: reconstruct
+ public :: init_reconstruct, reconstruct
contains
- subroutine reconstruct(s, grid)
+ subroutine init_reconstruct(grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: reconstruct vector field at cell centers based on radial basis functions
!
@@ -19,14 +19,8 @@
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
+ type (grid_meta), intent(inout) :: grid
- ! 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
@@ -39,18 +33,21 @@
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
+ real (kind=RKIND), dimension(:,:,:), pointer :: matrix_reconstruct
+ real (kind=RKIND), dimension(:,:,:), pointer :: normal
+ real (kind=RKIND), dimension(:,:), pointer :: rbf_value
+ real (kind=RKIND), allocatable, dimension(:,:) :: mwork
!========================================================
- ! needed for init procedure
+ ! arrays filled and saved during init procedure
!========================================================
+ matrix_reconstruct => grid % matrix_reconstruct % array
+ normal => grid % normal % array
+ rbf_value => grid % rbf_value % array
+
+ !========================================================
+ ! temporary variables needed for init procedure
+ !========================================================
xCell => grid % xCell % array
yCell => grid % yCell % array
zCell => grid % zCell % array
@@ -64,31 +61,22 @@
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
+ ! allocate arrays
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))
+
+ ! 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
- ! 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
+ ! 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)
@@ -107,64 +95,108 @@
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
+ 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)
- ! 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
+ enddo ! iCell
- ! 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)
+ deallocate(ipvt)
+ deallocate(work)
+ deallocate(xLoc)
- enddo ! iCell
+ end subroutine init_reconstruct
-!=================================================================
-! this ends the init procedure, the compute procedure begins below
-!=================================================================
+ 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
+ integer :: iCell,iEdge, i, j, k, npoints, EdgeMax
+ real (kind=RKIND) :: r
+ 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
+
+ ! stored arrays used during compute procedure
+ matrix_reconstruct => grid % matrix_reconstruct % array
+ normal => grid % normal % array
+ rbf_value => grid % rbf_value % array
+
+ ! temporary variables
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+ 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
@@ -177,12 +209,11 @@
! 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
+ rhs = 0.0
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
+ rhs(k,j) = u(k,iEdge)
enddo
enddo
@@ -196,7 +227,7 @@
enddo
enddo
- ! recontruct the velocity field at point X1 (Eq 6 in Tex document)
+ ! 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)
@@ -205,29 +236,13 @@
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
deallocate(rhs)
deallocate(c)
- ! deallocate (these deallocates will be removed when init procedure is implemented)
- deallocate(normal)
- deallocate(matrix_reconstruct)
- deallocate(rbf_value)
-
- contains
+ end subroutine reconstruct
subroutine get_distance(x1,x2,xout)
implicit none
@@ -261,16 +276,11 @@
! xout = exp(-r**2)
! multiquadrics
- ! xout = 1.0 / sqrt(1.0**2 + r**2)
+ xout = 1.0 / sqrt(1.0**2 + r**2)
! other
- xout = 1.0 / (1.0**2 + r**2)
+ ! xout = 1.0 / (1.0**2 + r**2)
end subroutine evaluate_rbf
-
-
- end subroutine reconstruct
-
-
end module vector_reconstruction
</font>
</pre>