<p><b>xylar@lanl.gov</b> 2010-04-27 12:41:17 -0600 (Tue, 27 Apr 2010)</p><p>Branch Commit<br>
<br>
More interpolation test subroutines -- this time for scalar RBFs with<br>
constant or constant+linear function exact reconstruction in 3D and in<br>
a 3D plane (2D in a unique coordinate system for each cell)<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-04-27 15:43:25 UTC (rev 214)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-04-27 18:41:17 UTC (rev 215)
@@ -26,9 +26,12 @@
print *, "2D RBF test"
call test_ibInterp_get2DRBF
-print *, "ScalarRBFAndConstant test"
- call test_ibInterp_computeScalarRBFAndConstant
+print *, "ScalarRBF test"
+ call test_ibInterp_computeScalarRBF
+print *, "ScalarRBFAndPlanar test"
+ call test_ibInterp_computeScalarRBFAndPlane
+
print *, "tests finished, stopping..."
stop
@@ -125,7 +128,7 @@
end subroutine test_ibInterp_get2DRBF
- subroutine test_ibInterp_computeScalarRBFAndConstant
+ subroutine test_ibInterp_computeScalarRBF
integer, parameter :: pointCount = 32
real(kind=RKIND), dimension(pointCount,3) :: sourcePoints
@@ -140,6 +143,8 @@
real(kind=RKIND) :: epsilon, alpha
+ integer :: testType
+
epsilon = 1e-6
alpha = 1.35
@@ -189,10 +194,18 @@
*sourcePoints(27:32,2) + interfaceNormals(27:32,2)*sourcePoints(27:32,1)**2
print *, "testFunction:", testFunction
+do testType = 1,2
+print *, "testType: ", testType
destinationPoint = sourcePoints(13,:)
- call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
- sourcePoints, isInterface, interfaceNormals, destinationPoint, &
- alpha, dirichletCoefficients, neumannCoefficients)
+ if(testType == 1) then
+ call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletCoefficients, neumannCoefficients)
+ else if(testType == 2) then
+ call ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletCoefficients, neumannCoefficients)
+ end if
functionAtDestination(1) = sum(testFunction*neumannCoefficients)
@@ -201,9 +214,15 @@
destinationPoint = sourcePoints(32,:) + epsilon*interfaceNormals(1,:)
- call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
- sourcePoints, isInterface, interfaceNormals, destinationPoint, &
- alpha, dirichletCoefficients, neumannCoefficients)
+ if(testType == 1) then
+ call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletCoefficients, neumannCoefficients)
+ else if(testType == 2) then
+ call ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletCoefficients, neumannCoefficients)
+ end if
!print *, "dirichlet:", dirichletCoefficients
!print *, "neumann:", neumannCoefficients
@@ -214,175 +233,115 @@
destinationPoint = sourcePoints(32,:) - epsilon*interfaceNormals(1,:)
- call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
- sourcePoints, isInterface, interfaceNormals, destinationPoint, &
- alpha, dirichletCoefficients, neumannCoefficients)
+ if(testType == 1) then
+ call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletCoefficients, neumannCoefficients)
+ else if(testType == 2) then
+ call ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletCoefficients, neumannCoefficients)
+ end if
functionAtDestination(2) = sum(testFunction*neumannCoefficients)
print *, "function value:", functionAtDestination
print *, "function normal deriv:", (functionAtDestination(1)-functionAtDestination(2))/(2.0*epsilon)
print *, "deriv should be:", testFunction(32)
+end do
- end subroutine test_ibInterp_computeScalarRBFAndConstant
+ end subroutine test_ibInterp_computeScalarRBF
-! subroutine test_ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
-! sourcePoints, isInterface, interfaceNormals, destinationPoint, &
-! alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
+ subroutine test_ibInterp_computeScalarRBFAndPlane
+ integer, parameter :: pointCount = 12
+ real(kind=RKIND), dimension(pointCount,3) :: sourcePoints
+ logical, dimension(pointCount) :: isInterface
+ real(kind=RKIND), dimension(pointCount,3) :: interfaceNormals
+ real(kind=RKIND), dimension(3) :: destinationPoint
+ real(kind=RKIND), dimension(pointCount) :: &
+ dirichletCoefficients, neumannCoefficients
+
+ real(kind=RKIND), dimension(pointCount) :: testFunction
+ real(kind=RKIND), dimension(2) :: functionAtDestination
+ real(kind=RKIND), dimension(2,3) :: planarBasisVectors
+
+ real(kind=RKIND) :: epsilon, alpha
+
+ epsilon = 1e-6
+ alpha = 1.35
+
+ ! set up a fake set of fluid points, interface points and interfce normals
+ sourcePoints(1,:) = (/0.0, -0.707106781186548, 0.707106781186548/)
+ sourcePoints(2,:) = (/0.816496580927726, -1.11535507165041, 0.298858490722685/)
+ sourcePoints(3,:) = (/-0.816496580927726, 0.408248290463863, 0.408248290463863/)
+ sourcePoints(4,:) = (/0.0, 0.0, 0.0/)
+ sourcePoints(5,:) = (/0.816496580927726, -0.408248290463863, -0.408248290463863/)
+ sourcePoints(6,:) = (/1.63299316185545, -0.816496580927726, -0.816496580927726/)
+ sourcePoints(7,:) = (/-0.816496580927726, 1.11535507165041, -0.298858490722685/)
+ sourcePoints(8,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
+ sourcePoints(9,:) = (/0.816496580927726, 0.298858490722685, -1.11535507165041/)
+ sourcePoints(10,:) = (/1.46969384566991, -0.310582854123025, -1.15911099154688/)
+ sourcePoints(11,:) = (/-0.163299316185545, 1.21302050799125, -1.04972119180570/)
+ sourcePoints(12,:) = (/0.489897948556636, 0.603579163145540, -1.09347711170218/)
+
+ interfaceNormals(:,1) = 0.365148371670111
+ interfaceNormals(:,2) = 0.449881346198621
+ interfaceNormals(:,3) = -0.815029717868732
+
+ isInterface(1:9) = .false.
+ isInterface(10:12) = .true.
+
+ testFunction(1) = 0
+ testFunction(2) = -1
+ testFunction(3) = 0
+ testFunction(4) = 0
+ testFunction(5) = 0
+ testFunction(6) = 0
+ testFunction(7) = 1
+ testFunction(8) = 0
+ testFunction(9) = 1
+ testFunction(10) = 3.86392546511964
+ testFunction(11) = -0.250439613479976
+ testFunction(12) = 0.965981366279909
+
+
+ planarBasisVectors(1,:) = (/0.816496580927726, -0.408248290463863, -0.408248290463863/)
+ planarBasisVectors(2,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
+
+ destinationPoint = sourcePoints(9,:)
+ call ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, planarBasisVectors, dirichletCoefficients, neumannCoefficients)
+
+ functionAtDestination(1) = sum(testFunction*neumannCoefficients)
+
+print *, "function value:", functionAtDestination(1)
+print *, "function should be:", testFunction(9)
+
+ destinationPoint = sourcePoints(12,:) + epsilon*interfaceNormals(1,:)
+
+ call ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, planarBasisVectors, dirichletCoefficients, neumannCoefficients)
+
+!print *, "dirichlet:", dirichletCoefficients
+!print *, "neumann:", neumannCoefficients
+
+ ! test to make sure it reconstructs a function properly
+
+ functionAtDestination(1) = sum(testFunction*neumannCoefficients)
+
+ destinationPoint = sourcePoints(12,:) - epsilon*interfaceNormals(1,:)
+
+ call ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, planarBasisVectors, dirichletCoefficients, neumannCoefficients)
+
+ functionAtDestination(2) = sum(testFunction*neumannCoefficients)
+print *, "function value:", functionAtDestination
+print *, "function normal deriv:", (functionAtDestination(1)-functionAtDestination(2))/(2.0*epsilon)
+print *, "deriv should be:", testFunction(12)
+ end subroutine test_ibInterp_computeScalarRBFAndPlane
!
-! integer, intent(in) :: pointCount
-! real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
-! logical, dimension(pointCount), intent(in) :: isInterface
-! real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
-! real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
-! real(kind=RKIND), intent(in) :: alpha
-! real(kind=RKIND), dimension(2,3) :: planeBasisVectors
-! real(kind=RKIND), dimension(pointCount), intent(out) :: &
-! dirichletCoefficients, neumannCoefficients
-! type (dm_info), intent(in) :: dminfo
-!
-! integer :: i, j
-! integer :: matrixSize
-!
-! real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
-! real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
-! integer, dimension(:), pointer :: pivotIndices
-!
-! matrixSize = pointCount+3 !! 3 extra space for constant and 2 planar dimensions
-!
-! allocate(dirichletMatrix(matrixSize,matrixSize))
-! allocate(neumannMatrix(matrixSize,matrixSize))
-! allocate(rhs(matrixSize))
-! allocate(rhsCopy(matrixSize))
-! allocate(coeffs(matrixSize))
-! allocate(pivotIndices(matrixSize))
-!
-! dirichletMatrix = 0.0
-! neumannMatrix = 0.0
-! rhs = 0.0
-! rhsCopy = 0.0
-! coeffs = 0.0
-!
-! call setUpScalarRBFMatrixAndRHS(pointCount, &
-! sourcePoints, isInterface, interfaceNormals, destinationPoint, &
-! alpha, dirichletMatrix(1:pointCount,1:pointCount), &
-! neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
-!
-! do i = 1, pointCount
-! dirichletMatrix(i,pointCount+1) = 1.0
-! dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVectors(1,:))
-! dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVectors(2,:))
-! if(isInterface(i)) then
-! neumannMatrix(i,pointCount+1) = 0.0
-! neumannMatrix(i,pointCount+2) = sum(interfaceNormals(i,1:3)*planeBasisVectors(1,:))
-! neumannMatrix(i,pointCount+3) = sum(interfaceNormals(i,1:3)*planeBasisVectors(2,:))
-! else
-! neumannMatrix(i,pointCount+1:pointCount+3) &
-! = dirichletMatrix(i,pointCount+1:pointCount+3)
-! end if
-! dirichletMatrix(pointCount+1:pointCount+3,i) &
-! = dirichletMatrix(i,pointCount+1:pointCount+3)
-! neumannMatrix(pointCount+1:pointCount+3,i) &
-! = neumannMatrix(i,pointCount+1:pointCount+3)
-! end do
-!
-! rhs(pointCount+1) = 1.0
-! rhs(pointCount+2) = sum(destinationPoint(1:3)*planeBasisVectors(1,:))
-! rhs(pointCount+3) = sum(destinationPoint(1:3)*planeBasisVectors(2,:))
-!
-! ! solve each linear system
-! rhsCopy = rhs
-! call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
-! dirichletCoefficients = coeffs(1:pointCount)
-!
-! call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
-! neumannCoefficients = coeffs(1:pointCount)
-!
-! deallocate(dirichletMatrix)
-! deallocate(neumannMatrix)
-! deallocate(rhs)
-! deallocate(rhsCopy)
-! deallocate(coeffs)
-! deallocate(pivotIndices)
-!
-! end subroutine test_ibInterp_computeScalarRBFAndPlaneCoefficients
-!
-! subroutine test_ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
-! sourcePoints, isInterface, interfaceNormals, destinationPoint, &
-! alpha, dirichletCoefficients, neumannCoefficients, dminfo)
-!
-! integer, intent(in) :: pointCount
-! real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
-! logical, dimension(pointCount), intent(in) :: isInterface
-! real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
-! real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
-! real(kind=RKIND), intent(in) :: alpha
-! real(kind=RKIND), dimension(pointCount), intent(out) :: &
-! dirichletCoefficients, neumannCoefficients
-! type (dm_info), intent(in) :: dminfo
-!
-! integer :: i, j
-! integer :: matrixSize
-!
-! real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
-! real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
-! integer, dimension(:), pointer :: pivotIndices
-!
-! matrixSize = pointCount+4 !! 4 extra space for constant and linear in 3D
-!
-! allocate(dirichletMatrix(matrixSize,matrixSize))
-! allocate(neumannMatrix(matrixSize,matrixSize))
-! allocate(rhs(matrixSize))
-! allocate(rhsCopy(matrixSize))
-! allocate(coeffs(matrixSize))
-! allocate(pivotIndices(matrixSize))
-!
-! dirichletMatrix = 0.0
-! neumannMatrix = 0.0
-! rhs = 0.0
-! rhsCopy = 0.0
-! coeffs = 0.0
-!
-! call setUpScalarRBFMatrixAndRHS(pointCount, &
-! sourcePoints, isInterface, interfaceNormals, destinationPoint, &
-! alpha, dirichletMatrix(1:pointCount,1:pointCount), &
-! neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
-!
-! do i = 1, pointCount
-! dirichletMatrix(i,pointCount+1) = 1.0
-! dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
-! if(isInterface(i)) then
-! neumannMatrix(i,pointCount+1) = 0.0
-! neumannMatrix(i,pointCount+2:pointCount+4) = interfaceNormals(i,1:3)
-! else
-! neumannMatrix(i,pointCount+1:pointCount+4) &
-! = dirichletMatrix(i,pointCount+1:pointCount+4)
-! end if
-! dirichletMatrix(pointCount+1:pointCount+4,i) &
-! = dirichletMatrix(i,pointCount+1:pointCount+4)
-! neumannMatrix(pointCount+1:pointCount+4,i) &
-! = neumannMatrix(i,pointCount+1:pointCount+4)
-! end do
-!
-! rhs(pointCount+1) = 1.0
-! rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
-!
-! ! solve each linear system
-! rhsCopy = rhs
-! call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
-! dirichletCoefficients = coeffs(1:pointCount)
-!
-! call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
-! neumannCoefficients = coeffs(1:pointCount)
-!
-! deallocate(dirichletMatrix)
-! deallocate(neumannMatrix)
-! deallocate(rhs)
-! deallocate(rhsCopy)
-! deallocate(coeffs)
-! deallocate(pivotIndices)
-!
-! end subroutine test_ibInterp_computeScalarRBFAndLinearCoefficients
-!
! subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &
! sourcePoints, unitVectors, destinationPoint, &
! alpha, coefficients, dminfo)
</font>
</pre>