<p><b>duda</b> 2010-03-30 14:01:26 -0600 (Tue, 30 Mar 2010)</p><p>Add vector reconstruction module as a shared (i.e., common to all<br>
cores) module in a new directory, src/operators. Make necessary<br>
changes to Makefiles to accommodate this new directory, and remove<br>
references to BLAS and LAPACK, which are no longer needed by any<br>
code.<br>
<br>
<br>
M src/core_hyd_atmos/Registry<br>
M src/core_hyd_atmos/Makefile<br>
M src/core_sw/mpas_interface.F<br>
M src/core_sw/Registry<br>
D src/core_sw/module_vector_reconstruction.F<br>
M src/core_sw/Makefile<br>
M src/core_ocean/Registry<br>
D src/core_ocean/module_vector_reconstruction.F<br>
M src/core_ocean/Makefile<br>
M src/Makefile<br>
A src/operators<br>
A src/operators/module_vector_reconstruction.F<br>
A src/operators/Makefile<br>
M Makefile<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/Makefile
===================================================================
--- trunk/mpas/Makefile 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/Makefile 2010-03-30 20:01:26 UTC (rev 171)
@@ -85,7 +85,7 @@
CPPINCLUDES = -I../inc -I$(NETCDF)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include
-LIBS = -L$(NETCDF)/lib -lnetcdf -L$(BLAS)/lib -L$(LAPACK)/lib -lblas -llapack
+LIBS = -L$(NETCDF)/lib -lnetcdf
RM = rm -f
CPP = cpp -C -P -traditional
Modified: trunk/mpas/src/Makefile
===================================================================
--- trunk/mpas/src/Makefile 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/Makefile 2010-03-30 20:01:26 UTC (rev 171)
@@ -2,8 +2,8 @@
all: mpas
-mpas: reg_includes externals frame dycore drver
- $(FC) $(LDFLAGS) -o $(CORE)_model.exe driver/*.o -L. -ldycore -lframework $(LIBS)
+mpas: reg_includes externals frame ops dycore drver
+ $(FC) $(LDFLAGS) -o $(CORE)_model.exe driver/*.o -L. -ldycore -lops -lframework $(LIBS)
reg_includes:
( cd registry; make CC="$(SCC)" )
@@ -16,6 +16,10 @@
( cd framework; make all )
ln -sf framework/libframework.a libframework.a
+ops:
+ ( cd operators; make all )
+ ln -sf operators/libops.a libops.a
+
dycore:
( cd core_$(CORE); make all )
ln -sf core_$(CORE)/libdycore.a libdycore.a
@@ -24,10 +28,11 @@
( cd driver; make all )
clean:
- $(RM) $(CORE)_model.exe libframework.a libdycore.a
+ $(RM) $(CORE)_model.exe libframework.a libops.a libdycore.a
( cd registry; make clean )
( cd external; make clean )
( cd framework; make clean )
+ ( cd operators; make clean )
( cd inc; rm -f *.inc )
if [ -d core_$(CORE) ] ; then \
( cd core_$(CORE); make clean ) \
Modified: trunk/mpas/src/core_hyd_atmos/Makefile
===================================================================
--- trunk/mpas/src/core_hyd_atmos/Makefile 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_hyd_atmos/Makefile 2010-03-30 20:01:26 UTC (rev 171)
@@ -24,5 +24,5 @@
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
- $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
+ $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
$(RM) $*.f90
Modified: trunk/mpas/src/core_hyd_atmos/Registry
===================================================================
--- trunk/mpas/src/core_hyd_atmos/Registry 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_hyd_atmos/Registry 2010-03-30 20:01:26 UTC (rev 171)
@@ -36,6 +36,7 @@
dim vertexDegree vertexDegree
dim FIFTEEN 15
dim TWENTYONE 21
+dim R3 3
dim nVertLevels nVertLevels
#dim nTracers nTracers
dim nVertLevelsP1 nVertLevels+1
@@ -123,6 +124,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 - -
@@ -152,3 +156,6 @@
var real deriv_two ( FIFTEEN TWO nEdges ) o deriv_two - -
var integer advCells ( TWENTYONE nCells ) - advCells - -
+# Arrays required for reconstruction of velocity field
+var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
+
Modified: trunk/mpas/src/core_ocean/Makefile
===================================================================
--- trunk/mpas/src/core_ocean/Makefile 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_ocean/Makefile 2010-03-30 20:01:26 UTC (rev 171)
@@ -2,7 +2,6 @@
OBJS = module_test_cases.o \
module_time_integration.o \
- module_vector_reconstruction.o \
module_global_diagnostics.o \
mpas_interface.o
@@ -13,13 +12,11 @@
module_test_cases.o:
-module_time_integration.o: module_vector_reconstruction.o
+module_time_integration.o:
-module_vector_reconstruction.o:
-
module_global_diagnostics.o:
-mpas_interface.o: module_global_diagnostics.o module_vector_reconstruction.o module_test_cases.o module_time_integration.o
+mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o
clean:
$(RM) *.o *.mod libdycore.a
@@ -27,5 +24,5 @@
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
- $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
+ $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
$(RM) $*.f90
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_ocean/Registry 2010-03-30 20:01:26 UTC (rev 171)
@@ -82,9 +82,6 @@
var real rho ( nVertLevels nCells Time ) iro rho - -
# 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
Deleted: trunk/mpas/src/core_ocean/module_vector_reconstruction.F
===================================================================
--- trunk/mpas/src/core_ocean/module_vector_reconstruction.F 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_ocean/module_vector_reconstruction.F 2010-03-30 20:01:26 UTC (rev 171)
@@ -1,354 +0,0 @@
-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
- ! xsad 10-02-09:
- !!!!! the lapack stuff isn't working on coyote in parallel !!!!!!!
- !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)
- ! xsad 10-02-09 end
-
- ! 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
-
- 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)
- 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
-!
- ! a more efficient reconstruction where rbf_values*matrix_reconstruct has been precomputed
- ! in coeffs_reconstruct
- 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
-
-
- enddo ! iCell
-
- ! deallocate
- deallocate(rhs)
-! deallocate(c)
-
- 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
Modified: trunk/mpas/src/core_sw/Makefile
===================================================================
--- trunk/mpas/src/core_sw/Makefile 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_sw/Makefile 2010-03-30 20:01:26 UTC (rev 171)
@@ -1,7 +1,6 @@
.SUFFIXES: .F .o
OBJS = module_test_cases.o \
- module_vector_reconstruction.o \
module_time_integration.o \
mpas_interface.o
@@ -12,17 +11,15 @@
module_test_cases.o:
-module_vector_reconstruction.o:
+module_time_integration.o:
-module_time_integration.o: module_vector_reconstruction.o
+mpas_interface.o: module_test_cases.o module_time_integration.o
-mpas_interface.o: module_test_cases.o module_time_integration.o module_vector_reconstruction.o
-
clean:
$(RM) *.o *.mod libdycore.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
- $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
+ $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
$(RM) $*.f90
Modified: trunk/mpas/src/core_sw/Registry
===================================================================
--- trunk/mpas/src/core_sw/Registry 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_sw/Registry 2010-03-30 20:01:26 UTC (rev 171)
@@ -78,9 +78,7 @@
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 - -
Deleted: trunk/mpas/src/core_sw/module_vector_reconstruction.F
===================================================================
--- trunk/mpas/src/core_sw/module_vector_reconstruction.F 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_sw/module_vector_reconstruction.F 2010-03-30 20:01:26 UTC (rev 171)
@@ -1,286 +0,0 @@
-module vector_reconstruction
-
- use grid_types
- use configure
- use constants
-
- 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), 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
-
- !========================================================
- ! 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
-
- ! 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,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
-
- 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
-
- 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
- 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
- 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
- rhs = 0.0
- do j=1,npoints
- iEdge = edgesOnCell(j,iCell)
- do k=1,nVertLevels
- rhs(k,j) = u(k,iEdge)
- 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
-
- enddo ! iCell
-
- ! deallocate
- deallocate(rhs)
- deallocate(c)
-
- 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(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 module vector_reconstruction
Modified: trunk/mpas/src/core_sw/mpas_interface.F
===================================================================
--- trunk/mpas/src/core_sw/mpas_interface.F 2010-03-30 18:06:43 UTC (rev 170)
+++ trunk/mpas/src/core_sw/mpas_interface.F 2010-03-30 20:01:26 UTC (rev 171)
@@ -26,6 +26,7 @@
call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh)
call init_reconstruct(mesh)
+ call reconstruct(block % time_levs(1) % state, mesh)
end subroutine mpas_init
Added: trunk/mpas/src/operators/Makefile
===================================================================
--- trunk/mpas/src/operators/Makefile (rev 0)
+++ trunk/mpas/src/operators/Makefile 2010-03-30 20:01:26 UTC (rev 171)
@@ -0,0 +1,19 @@
+.SUFFIXES: .F .o
+
+OBJS = module_vector_reconstruction.o
+
+all: operators
+
+operators: $(OBJS)
+ ar -ru libops.a $(OBJS)
+
+module_vector_reconstruction.o:
+
+clean:
+ $(RM) *.o *.mod libops.a
+
+.F.o:
+ $(RM) $@ $*.mod
+ $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
+ $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
+ $(RM) $*.f90
Added: trunk/mpas/src/operators/module_vector_reconstruction.F
===================================================================
--- trunk/mpas/src/operators/module_vector_reconstruction.F (rev 0)
+++ trunk/mpas/src/operators/module_vector_reconstruction.F 2010-03-30 20:01:26 UTC (rev 171)
@@ -0,0 +1,454 @@
+module vector_reconstruction
+
+ use grid_types
+ use configure
+ use constants
+
+ implicit none
+
+ public :: init_reconstruct, reconstruct
+
+ contains
+
+ subroutine init_reconstruct(grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: pre-compute coefficients used by the reconstruct() routine
+ !
+ ! Input: grid meta data
+ !
+ ! Output: grid % coeffs_reconstruct - coefficients used to reconstruct
+ ! velocity vectors 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, matrixSize
+ integer :: lwork, info
+ integer, allocatable, dimension(:) :: pivotingIndices
+ real (kind=RKIND), dimension(:), pointer :: dcEdge, xCell, yCell, zCell
+ real (kind=RKIND) :: r, rbfValue, v, X1(3), X2(3), alpha, rHat(3), normalDotRHat
+ real (kind=RKIND) :: xPlane, yPlane, xNormalPlane, yNormalPlane, xHatPlane(3), yHatPlane(3)
+ real (kind=RKIND), allocatable, dimension(:,:,:) :: xLoc
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+ real (kind=RKIND), allocatable, dimension(:,:) :: mwork
+ real (kind=RKIND), dimension(:,:), pointer :: matrix, invMatrix
+ real (kind=RKIND), dimension(:,:), pointer :: normals
+
+ !========================================================
+ ! arrays filled and saved during init procedure
+ !========================================================
+ 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(normals(3,EdgeMax))
+
+ ! init arrays
+ coeffs_reconstruct = 0.0
+ normals = 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
+ cell1 = iCell
+ X1(1) = xCell(cell1)
+ X1(2) = yCell(cell1)
+ X1(3) = zCell(cell1)
+
+ rHat = X1
+ call unit_vector_in_R3(rHat)
+
+ do j=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(j,iCell)
+ if (iCell == cellsOnEdge(1,iEdge)) then
+ cell2 = cellsOnEdge(2,iEdge)
+ X2(1) = xCell(cell2)
+ X2(2) = yCell(cell2)
+ X2(3) = zCell(cell2)
+ normals(:,j) = X2(:) - X1(:)
+ else
+ cell2 = cellsOnEdge(1,iEdge)
+ X2(1) = xCell(cell2)
+ X2(2) = yCell(cell2)
+ X2(3) = zCell(cell2)
+ normals(:,j) = X1(:) - X2(:)
+ endif
+
+ call unit_vector_in_R3(normals(:,j))
+
+ xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
+ enddo
+
+ npoints = nEdgesOnCell(iCell) ! only loop over the number of edges for this cell
+ matrixSize = npoints+2 ! room for 2 vector components for constant flow
+ ! basis functions
+ allocate(matrix(matrixSize,matrixSize))
+ matrix = 0.0
+ 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,rbfValue)
+ call get_dotproduct(normals(:,i),normals(:,j),v)
+ matrix(i,j) = v*rbfValue
+ enddo
+ enddo
+
+ ! add matrix entries to suppoert constant flow
+ ! xHat and yHat are a local basis in the plane of the horizontal cell
+ ! we arbitrarily choose xHat to point toward the first edge
+ call get_dotproduct(normals(:,1),rHat,normalDotRHat)
+ xHatPlane = normals(:,1) - normalDotRHat*rHat(:)
+ call unit_vector_in_R3(xHatPlane)
+
+ call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
+ call unit_vector_in_R3(yHatPlane) ! just to be sure...
+
+ do j=1,npoints
+ call get_dotproduct(normals(:,j),xHatPlane, xNormalPlane)
+ call get_dotproduct(normals(:,j),yHatPlane, yNormalPlane)
+ matrix(j,npoints+1) = xNormalPlane
+ matrix(j,npoints+2) = yNormalPlane
+ matrix(npoints+1,j) = matrix(j,npoints+1)
+ matrix(npoints+2,j) = matrix(j,npoints+2)
+ end do
+
+
+ ! invert matrix
+ allocate(invMatrix(matrixSize,matrixSize))
+ allocate(pivotingIndices(matrixSize))
+ invMatrix = 0.0
+ pivotingIndices = 0
+ call migs(matrix,matrixSize,invMatrix,pivotingIndices)
+
+ ! 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 i=1,npoints
+ ! compute value of RBF when evaluated between reconstruction location and edge locations
+ call get_distance(xLoc(:,i,iCell), X1(:), r)
+ r = r / alpha
+ call evaluate_rbf(r,rbfValue)
+
+ ! project the normals onto tangent plane of the cell
+ ! normal = normal - (normal dot r_hat) r_hat
+ ! rHat, the unit vector pointing from the domain center to the cell center
+ call get_dotproduct(normals(:,i),rHat,normalDotRHat)
+ normals(:,i) = normals(:,i) - normalDotRHat*rHat(:)
+
+ do j=1,npoints
+ coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
+ + rbfValue * normals(:,i) * invMatrix(i,j)
+ enddo
+ enddo
+ ! polynomial parts
+ do j=1,npoints
+ coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
+ + invMatrix(npoints+1,j)*xHatPlane
+ coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
+ + invMatrix(npoints+2,j)*yHatPlane
+ enddo
+
+ deallocate(matrix)
+ deallocate(invMatrix)
+ deallocate(pivotingIndices)
+
+ enddo ! iCell
+
+ deallocate(xLoc)
+ deallocate(normals)
+
+ 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 :: nCellsSolve
+ integer, dimension(:,:), pointer :: edgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer :: iCell,iEdge, i
+ real (kind=RKIND), dimension(:,:), pointer :: u
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+
+ ! stored arrays used during compute procedure
+ coeffs_reconstruct => grid % coeffs_reconstruct % array
+
+ ! temporary variables
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ nCellsSolve = grid % nCellsSolve
+ u => s % u % array
+ uReconstructX => s % uReconstructX % array
+ uReconstructY => s % uReconstructY % array
+ uReconstructZ => s % uReconstructZ % array
+
+ ! init the intent(out)
+ uReconstructX = 0.0
+ uReconstructY = 0.0
+ uReconstructZ = 0.0
+
+ ! loop over cell centers
+ do iCell=1,nCellsSolve
+ ! a more efficient reconstruction where rbf_values*matrix_reconstruct has been precomputed
+ ! in coeffs_reconstruct
+ do i=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i,iCell)
+ uReconstructX(:,iCell) = uReconstructX(:,iCell) &
+ + coeffs_reconstruct(1,i,iCell) * u(:,iEdge)
+ uReconstructY(:,iCell) = uReconstructY(:,iCell) &
+ + coeffs_reconstruct(2,i,iCell) * u(:,iEdge)
+ uReconstructZ(:,iCell) = uReconstructZ(:,iCell) &
+ + coeffs_reconstruct(3,i,iCell) * u(:,iEdge)
+ enddo
+ enddo ! iCell
+
+ 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(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 + xin**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
+!======================================================================
+
+! Updated 10/24/2001.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! !
+! Please Note: !
+! !
+! (1) This computer program is written by Tao Pang in conjunction with !
+! his book, "An Introduction to Computational Physics," published !
+! by Cambridge University Press in 1997. !
+! !
+! (2) No warranties, express or implied, are made for this program. !
+! !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+SUBROUTINE MIGS (A,N,X,INDX)
+!
+! Subroutine to invert matrix A(N,N) with the inverse stored
+! in X(N,N) in the output. Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ INTEGER, INTENT (IN) :: N
+ INTEGER :: I,J,K
+ INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
+ REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
+ REAL (kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
+ REAL (kind=RKIND), DIMENSION (N,N) :: B
+!
+ DO I = 1, N
+ DO J = 1, N
+ B(I,J) = 0.0
+ END DO
+ END DO
+ DO I = 1, N
+ B(I,I) = 1.0
+ END DO
+!
+ CALL ELGS (A,N,INDX)
+!
+ DO I = 1, N-1
+ DO J = I+1, N
+ DO K = 1, N
+ B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K)
+ END DO
+ END DO
+ END DO
+!
+ DO I = 1, N
+ X(N,I) = B(INDX(N),I)/A(INDX(N),N)
+ DO J = N-1, 1, -1
+ X(J,I) = B(INDX(J),I)
+ DO K = J+1, N
+ X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
+ END DO
+ X(J,I) = X(J,I)/A(INDX(J),J)
+ END DO
+ END DO
+END SUBROUTINE MIGS
+
+
+SUBROUTINE ELGS (A,N,INDX)
+!
+! Subroutine to perform the partial-pivoting Gaussian elimination.
+! A(N,N) is the original matrix in the input and transformed matrix
+! plus the pivoting element ratios below the diagonal in the output.
+! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ INTEGER, INTENT (IN) :: N
+ INTEGER :: I,J,K,ITMP
+ INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
+ REAL (kind=RKIND) :: C1,PI,PI1,PJ
+ REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
+ REAL (kind=RKIND), DIMENSION (N) :: C
+!
+! Initialize the index
+!
+ DO I = 1, N
+ INDX(I) = I
+ END DO
+!
+! Find the rescaling factors, one from each row
+!
+ DO I = 1, N
+ C1= 0.0
+ DO J = 1, N
+ !C1 = AMAX1(C1,ABS(A(I,J)))
+ C1 = MAX(C1,ABS(A(I,J)))
+ END DO
+ C(I) = C1
+ END DO
+!
+! Search the pivoting (largest) element from each column
+!
+ DO J = 1, N-1
+ PI1 = 0.0
+ DO I = J, N
+ PI = ABS(A(INDX(I),J))/C(INDX(I))
+ IF (PI.GT.PI1) THEN
+ PI1 = PI
+ K = I
+ ENDIF
+ END DO
+!
+! Interchange the rows via INDX(N) to record pivoting order
+!
+ ITMP = INDX(J)
+ INDX(J) = INDX(K)
+ INDX(K) = ITMP
+ DO I = J+1, N
+ PJ = A(INDX(I),J)/A(INDX(J),J)
+!
+! Record pivoting ratios below the diagonal
+!
+ A(INDX(I),J) = PJ
+!
+! Modify other elements accordingly
+!
+ DO K = J+1, N
+ A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
+ END DO
+ END DO
+ END DO
+!
+END SUBROUTINE ELGS
+
+end module vector_reconstruction
</font>
</pre>