<p><b>xylar@lanl.gov</b> 2010-04-22 16:23:02 -0600 (Thu, 22 Apr 2010)</p><p></p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/IBinterp/mpas/Makefile
===================================================================
--- branches/ocean_projects/IBinterp/mpas/Makefile 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/Makefile 2010-04-22 22:23:02 UTC (rev 203)
@@ -75,14 +75,37 @@
"CC = mpicc" \
"SFC = gfortran" \
"SCC = gcc" \
- "FFLAGS = -O3 -m64 -ffree-line-length-none" \
+ "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3 -m64" \
"CORE = $(CORE)" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+g95:
+ ( make all \
+ "FC = mpif90" \
+ "CC = mpicc" \
+ "SFC = g95" \
+ "SCC = gcc" \
+ "FFLAGS = -O3 -ffree-line-length-huge -r8" \
+ "CFLAGS = -O3" \
+ "LDFLAGS = -O3" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+g95-serial:
+ ( make all \
+ "FC = g95" \
+ "CC = gcc" \
+ "SFC = g95" \
+ "SCC = gcc" \
+ "FFLAGS = -O3 -ffree-line-length-huge -r8" \
+ "CFLAGS = -O3" \
+ "LDFLAGS = -O3" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+
CPPINCLUDES = -I../inc -I$(NETCDF)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include
LIBS = -L$(NETCDF)/lib -lnetcdf
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/Registry 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/Registry 2010-04-22 22:23:02 UTC (rev 203)
@@ -70,6 +70,11 @@
var real areaCell ( nCells ) iro areaCell - -
var real areaTriangle ( nVertices ) iro areaTriangle - -
+var real edgeNormalVectors ( R3 nEdges ) o edgeNormalVectors - -
+var real radialUnitVectors ( R3 nCells ) o radialUnitVectors - -
+var real edgeLocations ( R3 nEdges ) o edgeLocations - -
+var real cellTangentPlane ( R3 TWO nEdges ) o cellTangentPlane - -
+
var integer cellsOnCell ( maxEdges nCells ) iro cellsOnCell - -
var integer verticesOnCell ( maxEdges nCells ) iro verticesOnCell - -
var integer verticesOnEdge ( TWO nEdges ) iro verticesOnEdge - -
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_test_cases.F 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_test_cases.F 2010-04-22 22:23:02 UTC (rev 203)
@@ -552,7 +552,7 @@
real (kind=RKIND), intent(in) :: theta
AA = 0.5 * w * (2.0 * omega + w) * cos(theta)**2.0 + &
- 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*R**2*cos(theta)**-2.0)
+ 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*R**2*cos(theta)**(-2.0))
end function AA
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_time_integration.F 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_time_integration.F 2010-04-22 22:23:02 UTC (rev 203)
@@ -44,7 +44,7 @@
! mrp 100310 I added this to avoid running with NaNs
if (isNaN(sum(block % time_levs(2) % state % u % array))) then
print *, 'Stopping: NaN detected'
- call MPI_abort(MPI_COMM_WORLD,errorcode,ierr)
+ call dmpar_abort(domain%dminfo)
endif
! mrp 100310 end
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-04-22 22:23:02 UTC (rev 203)
@@ -16,7 +16,10 @@
use grid_types
use time_integration
+ use IB_interpolation
+ use vector_reconstruction
+
implicit none
type (block_type), intent(inout) :: block
@@ -25,10 +28,8 @@
call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh)
- ! xsad 10-02-09:
- ! commenting this out until we incorporate the necessary lapack routines into mpas
- ! call init_reconstruct(block_ptr % mesh)
- ! xsad 10-02-09 end
+ call ibInterp_initialize(mesh)
+ call init_reconstruct(mesh)
! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
! input arguement into mpas_init. Ask about that later. For now, there will be
@@ -37,7 +38,6 @@
! call timer_start("global diagnostics")
! call computeGlobalDiagnostics(domain % dminfo, block % time_levs(1) % state, mesh, 0, dt)
! call timer_stop("global diagnostics")
-! ! xsad 10-02-08 end
! call output_state_init(output_obj, domain, "OUTPUT")
! call write_output_frame(output_obj, domain)
Modified: branches/ocean_projects/IBinterp/mpas/src/framework/module_io_input.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/framework/module_io_input.F 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/src/framework/module_io_input.F 2010-04-22 22:23:02 UTC (rev 203)
@@ -993,10 +993,10 @@
integer, dimension(:), pointer :: super_int1d
integer, dimension(:,:), pointer :: super_int2d
- real :: super_real0d
- real, dimension(:), pointer :: super_real1d
- real, dimension(:,:), pointer :: super_real2d
- real, dimension(:,:,:), pointer :: super_real3d
+ real (kind=RKIND) :: super_real0d
+ real (kind=RKIND), dimension(:), pointer :: super_real1d
+ real (kind=RKIND), dimension(:,:), pointer :: super_real2d
+ real (kind=RKIND), dimension(:,:,:), pointer :: super_real3d
integer :: k
Modified: branches/ocean_projects/IBinterp/mpas/src/framework/module_io_output.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/framework/module_io_output.F 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/src/framework/module_io_output.F 2010-04-22 22:23:02 UTC (rev 203)
@@ -126,9 +126,9 @@
integer, dimension(:), pointer :: super_int1d
integer, dimension(:,:), pointer :: super_int2d
real :: super_real0d
- real, dimension(:), pointer :: super_real1d
- real, dimension(:,:), pointer :: super_real2d
- real, dimension(:,:,:), pointer :: super_real3d
+ real (kind=RKIND), dimension(:), pointer :: super_real1d
+ real (kind=RKIND), dimension(:,:), pointer :: super_real2d
+ real (kind=RKIND), dimension(:,:,:), pointer :: super_real3d
output_obj % time = itime
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/Makefile
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/Makefile 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/Makefile 2010-04-22 22:23:02 UTC (rev 203)
@@ -1,6 +1,6 @@
.SUFFIXES: .F .o
-OBJS = module_vector_reconstruction.o
+OBJS = module_IB_interpolation.o module_vector_reconstruction.o
all: operators
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F 2010-04-22 21:22:15 UTC (rev 202)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F 2010-04-22 22:23:02 UTC (rev 203)
@@ -1,106 +1,114 @@
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-module IB_interpoation
+module IB_interpolation
use dmpar
+ use grid_types
implicit none
private
save
- integer, parameter, public :: &
+
+ integer, public, parameter :: &
MPAS_i4 = selected_int_kind(6), &
MPAS_r8 = selected_real_kind(13)
public :: ibInterp_initialize, ibInterp_get2DRBFCoefficients, &
ibInterp_evaluate2DRBFWithDerivs, &
- ibInterp_computeScalarRBFCoefficients, &
- ibInterp_computeVectorNoSlipCoefficientsRBF, &
- ibInterp_computeVectorFreeSlipCoefficientsRBF
+ ibInterp_computeScalarRBFAndConstantCoefficients, &
+ ibInterp_computeScalarRBFAndPlaneCoefficients, &
+ ibInterp_computeScalarRBFAndLinearCoefficients, &
+ ibInterp_computeVectorNoSlipRBFCoefficients, &
+ ibInterp_computeVectorNoSlipPlanarRBFCoefficients, &
+ ibInterp_computeVectorFreeSlipRBFCoefficients, &
+ ibInterp_computeVectorFreeSlipPlanarRBFCoefficients
integer(MPAS_i4), public, parameter :: &
kMPAS_IBRBFPolynomialConstant = 0, &
kMPAS_IBRBFPolynomialLinear = 1
-
contains
- subroutine ibInterp_initialize(grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose:
- !
- ! Input:
- !
- ! Output:
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine ibInterp_initialize(grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose:
+ !
+ ! Input:
+ !
+ ! Output:
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- implicit none
+ implicit none
- type (grid_meta), intent(inout) :: grid
+ type (grid_meta), intent(inout) :: grid
- integer(MPAS_i4) :: nCells, nEdges
- integer(MPAS_i4), dimension(:,:), pointer :: cellsOnEdge
- integer(MPAS_i4) :: iEdge, iCell, cell1, cell2
- real(MPAS_r8), dimension(:), pointer :: xCell, yCell, zCell
- real(MPAS_r8), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &
- edgeLocations
- real(MPAS_r8), dimension(3) :: xHatPlane, yHatPlane
- real(MPAS_r8) :: normalDotRHat
+ integer(MPAS_i4) :: nCells, nEdges
+ integer(MPAS_i4), dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
+ integer(MPAS_i4) :: iEdge, iCell, cell1, cell2
+ real(MPAS_r8), dimension(:), pointer :: xCell, yCell, zCell
+ real(MPAS_r8), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &
+ edgeLocations
+ real(MPAS_r8), dimension(:,:,:), pointer :: cellTangentPlane
+ real(MPAS_r8), dimension(3) :: xHatPlane, yHatPlane, rHat
+ real(MPAS_r8) :: normalDotRHat
- xCell => grid % xCell % array
- yCell => grid % yCell % array
- zCell => grid % zCell % array
- cellsOnEdge => grid % cellsOnEdge % array
- nCells = grid % nCells
- nEdges = grid % nEdges
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+ nCells = grid % nCells
+ nEdges = grid % nEdges
- radialUnitVectors => grid % radialUnitVectors % array
- edgeNormalVectors => grid % edgeNormalVectors % array
- edgeLocations => grid % edgeLocations % array
- cellTangentPlane => grid % cellTangentPlane % array
+ radialUnitVectors => grid % radialUnitVectors % array
+ edgeNormalVectors => grid % edgeNormalVectors % array
+ edgeLocations => grid % edgeLocations % array
+ cellTangentPlane => grid % cellTangentPlane % array
- ! init arrays
- edgeNormalVectors = 0
- radialUnitVectors = 0
+ ! init arrays
+ edgeNormalVectors = 0
+ radialUnitVectors = 0
- ! loop over all cells to be solved on this block
- do iCell=1,nCells
- radialUnitVectors(iCell,1) = xCell(iCell)
- radialUnitVectors(iCell,2) = yCell(iCell)
- radialUnitVectors(iCell,3) = zCell(iCell)
- call unit_vector_in_R3(radialUnitVectors(iCell,:))
- end do
+ ! loop over all cells to be solved on this block
+ do iCell=1,nCells
+ radialUnitVectors(1,iCell) = xCell(iCell)
+ radialUnitVectors(2,iCell) = yCell(iCell)
+ radialUnitVectors(3,iCell) = zCell(iCell)
+ call unit_vector_in_R3(radialUnitVectors(:,iCell))
+ end do
- do iEdge = 1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- edgeNormalVectors(iEdge,1) = xCell(cell2) - xCell(cell1)
- edgeNormalVectors(iEdge,2) = yCell(cell2) - yCell(cell1)
- edgeNormalVectors(iEdge,3) = zCell(cell2) - zCell(cell1)
- edgeLocations(iEdge,1) = 0.5*(xCell(cell2) + xCell(cell1))
- edgeLocations(iEdge,2) = 0.5*(yCell(cell2) + yCell(cell1))
- edgeLocations(iEdge,3) = 0.5*(zCell(cell2) + zCell(cell1))
- call unit_vector_in_R3(edgeNormalVectors(iEdge,:))
- end do
+ do iEdge = 1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ edgeNormalVectors(1,iEdge) = xCell(cell2) - xCell(cell1)
+ edgeNormalVectors(2,iEdge) = yCell(cell2) - yCell(cell1)
+ edgeNormalVectors(3,iEdge) = zCell(cell2) - zCell(cell1)
+ call unit_vector_in_R3(edgeNormalVectors(:,iEdge))
+ edgeLocations(1,iEdge) = 0.5*(xCell(cell2) + xCell(cell1))
+ edgeLocations(2,iEdge) = 0.5*(yCell(cell2) + yCell(cell1))
+ edgeLocations(3,iEdge) = 0.5*(zCell(cell2) + zCell(cell1))
+ end do
- do iCell=1,nCells
- iEdge = edgesOnCell(1,iCell)
- ! xHat and yHat are a local basis in the plane of the horizontal cell
- ! we arbitrarily choose xHat to point toward the first edge
- normalDotRHat = sum(edgeNormalVectors(iEdge,:)*radialUnitVectors(iCell,:))
- xHatPlane = edgeNormalVectors(iEdge,:) - normalDotRHat*radialUnitVectors(iCell,1)
- call unit_vector_in_R3(xHatPlane)
+ do iCell=1,nCells
+ iEdge = edgesOnCell(1,iCell)
+ ! xHat and yHat are a local basis in the plane of the horizontal cell
+ ! we arbitrarily choose xHat to point toward the first edge
+ rHat = radialUnitVectors(:,iCell)
+ normalDotRHat = sum(edgeNormalVectors(:,iEdge)*rHat)
+ xHatPlane = edgeNormalVectors(:,iEdge) - 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...
- cellTangentPlane(iCell, 1, :) = xHatPlane
- cellTangentPlane(iCell, 2, :) = yHatPlane
- end do
+ call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
+ call unit_vector_in_R3(yHatPlane) ! just to be sure...
+ cellTangentPlane(:,1,iCell) = xHatPlane
+ cellTangentPlane(:,2,iCell) = yHatPlane
+ end do
- end subroutine ibInterp_initialize
+ end subroutine ibInterp_initialize
- subroutine ibInterp_getRBFCoefficients(pointCount, coeffCount, dimensions, points, fieldValues, &
- polynomialFlag, alpha, coefficients, dminfo)
+ subroutine ibInterp_get2DRBFCoefficients(pointCount, coeffCount, dimensions, points, &
+ fieldValues, polynomialFlag, alpha, coefficients, dminfo)
integer(MPAS_i4), intent(in) :: pointCount, coeffCount, dimensions
real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: points
@@ -208,9 +216,9 @@
subroutine ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
- polynomialFlag, alpha, dirichletCoefficients, neumannCoefficients, dminfo)
+ alpha, dirichletCoefficients, neumannCoefficients, dminfo)
- integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+ integer(MPAS_i4), intent(in) :: pointCount
real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
logical, dimension(pointCount), intent(in) :: isInterface
real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
@@ -221,15 +229,24 @@
type (dm_info), intent(in) :: dminfo
integer(MPAS_i4) :: i, j
- integer(MPAS_i4), parameter :: matrixSize = pointCount+1 !! 1 extra space for constant
+ integer(MPAS_i4) :: matrixSize
+ real(MPAS_r8), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+ real(MPAS_r8), dimension(:), pointer :: rhs, rhsCopy, coeffs
+ integer(MPAS_i4), dimension(:), pointer :: pivotIndices
- real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
- real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
- integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+ matrixSize = pointCount+1 !! 1 extra space for constant
- real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+ allocate(dirichletMatrix(matrixSize,matrixSize))
+ allocate(neumannMatrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize))
+ allocate(rhsCopy(matrixSize))
+ allocate(coeffs(matrixSize))
+ allocate(pivotIndices(matrixSize))
+ dirichletMatrix = 0.0_MPAS_r8
+ neumannMatrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
rhsCopy = 0.0_MPAS_r8
coeffs = 0.0_MPAS_r8
@@ -259,13 +276,20 @@
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 ibInterp_computeScalarRBFAndConstantCoefficients
- subroutine ibInterp_computeScalarRBFAndPlane(pointCount, &
+ subroutine ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
- integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+ integer(MPAS_i4), intent(in) :: pointCount
real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
logical, dimension(pointCount), intent(in) :: isInterface
real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
@@ -277,14 +301,24 @@
type (dm_info), intent(in) :: dminfo
integer(MPAS_i4) :: i, j
- integer(MPAS_i4), parameter :: matrixSize = pointCount+3 ! extra space for constant and plane
+ integer(MPAS_i4) :: matrixSize
- real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
- real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
- integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+ real(MPAS_r8), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+ real(MPAS_r8), dimension(:), pointer :: rhs, rhsCopy, coeffs
+ integer(MPAS_i4), dimension(:), pointer :: pivotIndices
- real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+ 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_MPAS_r8
+ neumannMatrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
rhsCopy = 0.0_MPAS_r8
coeffs = 0.0_MPAS_r8
@@ -295,23 +329,25 @@
do i = 1, pointCount
dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
- dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVector(1,:))
- dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVector(2,:))
+ 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_MPAS_r8
- neumannMatrix(i,pointCount+2:pointCount+4) = interfaceNormals(i,1:3)
+ 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+4) &
- = dirichletMatrix(i,pointCount+1:pointCount+4)
+ neumannMatrix(i,pointCount+1:pointCount+3) &
+ = dirichletMatrix(i,pointCount+1:pointCount+3)
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)
+ 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_MPAS_r8
- rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
+ rhs(pointCount+2) = sum(destinationPoint(1:3)*planeBasisVectors(1,:))
+ rhs(pointCount+3) = sum(destinationPoint(1:3)*planeBasisVectors(2,:))
! solve each linear system
rhsCopy = rhs
@@ -321,32 +357,44 @@
call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
neumannCoefficients = coeffs(1:pointCount)
- end subroutine ibInterp_computeScalarRBFCoefficients
+ deallocate(dirichletMatrix)
+ deallocate(neumannMatrix)
+ deallocate(rhs)
+ deallocate(rhsCopy)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
- subroutine ibInterp_computeScalarRBFAndLinear(pointCount, &
+ end subroutine ibInterp_computeScalarRBFAndPlaneCoefficients
+
+ subroutine ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
- alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
+ alpha, dirichletCoefficients, neumannCoefficients, dminfo)
- integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+ integer(MPAS_i4), intent(in) :: pointCount
real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
logical, dimension(pointCount), intent(in) :: isInterface
real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
real(MPAS_r8), intent(in) :: alpha
- real(MPAS_r8), dimension(3,2) :: planeBasisVectors
real(MPAS_r8), dimension(pointCount), intent(out) :: &
dirichletCoefficients, neumannCoefficients
type (dm_info), intent(in) :: dminfo
integer(MPAS_i4) :: i, j
- integer(MPAS_i4), parameter :: matrixSize = pointCount+4 ! extra space for constant and plane
+ integer(MPAS_i4) :: matrixSize
+ real(MPAS_r8), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+ real(MPAS_r8), dimension(:), pointer :: rhs, rhsCopy, coeffs
+ integer(MPAS_i4), dimension(:), pointer :: pivotIndices
- real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
- real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
- integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+ matrixSize = pointCount+4 !! 4 extra space for constant and linear in 3D
- real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+ allocate(dirichletMatrix(matrixSize,matrixSize))
+ allocate(neumannMatrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize))
+ allocate(rhsCopy(matrixSize))
+ allocate(coeffs(matrixSize))
+ allocate(pivotIndices(matrixSize))
dirichletMatrix = 0.0_MPAS_r8
neumannMatrix = 0.0_MPAS_r8
@@ -354,70 +402,30 @@
rhsCopy = 0.0_MPAS_r8
coeffs = 0.0_MPAS_r8
- do j = 1, pointCount
- if(isInterface(j)) cycle
- do i = 1, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- rbfValue = evaluateRBF(rSquared)
- dirichletMatrix(i,j) = rbfValue
- neumannMatrix(i,j) = rbfValue
- end do
- end do
+ call setUpScalarRBFMatrixAndRHS(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletMatrix(1:pointCount,1:pointCount), &
+ neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
- do j = 1, pointCount
- if(.not. isInterface(j)) cycle
- do i = 1, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- normalDotX = sum(interfaceNormals(j,:) &
- * (sourcePoints(j,:)-sourcePoints(i,:)))
- call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
- rbfDerivOverR = rbfDerivOverR/alpha**2
- dirichletMatrix(i,j) = rbfValue
- neumannMatrix(i,j) = rbfDerivOverR*normalDotX
- end do
+ do i = 1, pointCount
+ dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+ dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
+ if(isInterface(i)) then
+ neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+ 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
- do j = 1, pointCount
- rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
- rhs(j) = evaluateRBF(rSquared)
- end do
+ rhs(pointCount+1) = 1.0_MPAS_r8
+ rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
- do i = 1, pointCount
- dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
- if(isInterface(i)) then
- neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
- else
- neumannMatrix(i,pointCount+1) = dirichletMatrix(i,pointCount+1)
- end if
- dirichletMatrix(pointCount+1,i) = dirichletMatrix(i,pointCount+1)
- neumannMatrix(pointCount+1,i) = neumannMatrix(i,pointCount+1)
- end do
-
- rhs(pointCount+1) = 1.0_MPAS_r8
- else if(polynomialFlag == kMPAS_IBRBFPolynomialLinear) then
- do i = 1, pointCount
- dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
- dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
- if(isInterface(i)) then
- neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
- 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_MPAS_r8
- rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
- else
- call exit_MPAS("ibInterp_computeScalarRBFCoefficients: not a valid polynomialFlag", dminfo)
- end if
-
-
! solve each linear system
rhsCopy = rhs
call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
@@ -426,194 +434,290 @@
call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
neumannCoefficients = coeffs(1:pointCount)
- end subroutine ibInterp_computeScalarRBFCoefficients
+ deallocate(dirichletMatrix)
+ deallocate(neumannMatrix)
+ deallocate(rhs)
+ deallocate(rhsCopy)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
- subroutine ibInterp_computeVectorNoSlipCoefficientsRBF(pointCount, &
- dimensions, sourcePoints, unitVectors, destinationPoint, &
+ end subroutine ibInterp_computeScalarRBFAndLinearCoefficients
+
+ subroutine ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &
+ sourcePoints, unitVectors, destinationPoint, &
alpha, coefficients, dminfo)
- integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
- real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
- real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
- real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
+ integer(MPAS_i4), intent(in) :: pointCount
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
real(MPAS_r8), intent(in) :: alpha
- real(MPAS_r8), dimension(pointCount, dimensions), intent(out) :: coefficients
+ real(MPAS_r8), dimension(pointCount, 3), intent(out) :: coefficients
type (dm_info), intent(in) :: dminfo
- integer(MPAS_i4) :: matrixSize, i, j
+ integer(MPAS_i4) :: i, j
+ integer(MPAS_i4) :: matrixSize
real(MPAS_r8), dimension(:,:), pointer :: matrix
real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
integer(MPAS_i4), dimension(:), pointer :: pivotIndices
- real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+ matrixSize = pointCount+3 ! extra space for constant vector
- matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,3))
+ allocate(coeffs(matrixSize,3))
+ allocate(pivotIndices(matrixSize))
- allocate(matrix(matrixSize,matrixSize))
- allocate(rhs(matrixSize,dimensions))
- allocate(coeffs(matrixSize,dimensions))
- allocate(pivotIndices(matrixSize))
-
matrix = 0.0_MPAS_r8
rhs = 0.0_MPAS_r8
coeffs = 0.0_MPAS_r8
- do j = 1, pointCount
- do i = 1, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
- rbfValue = evaluateRBF(rSquared)
- matrix(j,i) = rbfValue*unitVectorDotProduct
- end do
- rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
- rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
- end do
+ call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 3, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, matrix, rhs)
do i = 1, pointCount
- matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
- matrix(pointCount+1:pointCount+dimensions,i) &
- = matrix(i,pointCount+1:pointCount+dimensions)
+ matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
+ matrix(pointCount+1:pointCount+3,i) &
+ = matrix(i,pointCount+1:pointCount+3)
end do
- do i = 1, dimensions
- rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+ do i = 1, 3
+ rhs(pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
end do
! solve each linear system
- do i = 1, dimensions
+ do i = 1, 3
call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
end do
coefficients = coeffs(1:pointCount,:)
- end subroutine ibInterp_computeVectorNoSlipCoefficientsRBF
+ deallocate(matrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
- subroutine ibInterp_computeVectorNoSlipCoefficientsRBF(pointCount, &
- dimensions, sourcePoints, unitVectors, destinationPoint, &
- alpha, coefficients, dminfo)
+ end subroutine ibInterp_computeVectorNoSlipRBFCoefficients
- integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
- real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
- real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
- real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
+ subroutine ibInterp_computeVectorNoSlipPlanarRBFCoefficients(pointCount, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, planeBasisVectors, coefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
real(MPAS_r8), intent(in) :: alpha
- real(MPAS_r8), dimension(pointCount, dimensions), intent(out) :: coefficients
+ real(MPAS_r8), dimension(2,3) :: planeBasisVectors
+ real(MPAS_r8), dimension(pointCount, 3), intent(out) :: coefficients
type (dm_info), intent(in) :: dminfo
- integer(MPAS_i4) :: matrixSize, i, j
+ integer(MPAS_i4) :: i, j
+ integer(MPAS_i4) :: matrixSize
+ real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+
+ real(MPAS_r8), dimension(pointCount,2) :: planarSourcePoints
+ real(MPAS_r8), dimension(pointCount,2) :: planarUnitVectors
+ real(MPAS_r8), dimension(2) :: planarDestinationPoint
+
real(MPAS_r8), dimension(:,:), pointer :: matrix
real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
integer(MPAS_i4), dimension(:), pointer :: pivotIndices
- real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+ matrixSize = pointCount+2 ! space for constant vector in plane
- matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,2))
+ allocate(coeffs(matrixSize,2))
+ allocate(pivotIndices(matrixSize))
- allocate(matrix(matrixSize,matrixSize))
- allocate(rhs(matrixSize,dimensions))
- allocate(coeffs(matrixSize,dimensions))
- allocate(pivotIndices(matrixSize))
-
matrix = 0.0_MPAS_r8
rhs = 0.0_MPAS_r8
coeffs = 0.0_MPAS_r8
- do j = 1, pointCount
- do i = 1, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
- rbfValue = evaluateRBF(rSquared)
- matrix(j,i) = rbfValue*unitVectorDotProduct
- end do
- rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
- rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+ do i = 1, pointCount
+ planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:))
+ planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:))
+ planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+ planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
end do
+ planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:))
+ planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
+ call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 2, &
+ planarSourcePoints, planarUnitVectors, planarDestinationPoint, &
+ alpha, matrix, rhs)
do i = 1, pointCount
- matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
- matrix(pointCount+1:pointCount+dimensions,i) &
- = matrix(i,pointCount+1:pointCount+dimensions)
+ matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+ matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+ matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
end do
- do i = 1, dimensions
- rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+ rhs(pointCount+1,:) = planeBasisVectors(1,:)
+ rhs(pointCount+2,:) = planeBasisVectors(2,:)
+
+ ! solve each linear system
+ call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
+ call LEGS(matrix, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
+
+ coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,1)*coeffs(1:pointCount,2)
+ coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,2)*coeffs(1:pointCount,2)
+ coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,3)*coeffs(1:pointCount,2)
+
+ deallocate(matrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine ibInterp_computeVectorNoSlipPlanarRBFCoefficients
+
+ subroutine ibInterp_computeVectorFreeSlipRBFCoefficients(pointCount, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+ alpha, coefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isTangentToInterface
+ integer(MPAS_i4), dimension(pointCount), intent(in) :: normalVectorIndex
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(pointCount, 3), intent(out) :: coefficients
+ type (dm_info), intent(in) :: dminfo
+
+ integer(MPAS_i4) :: i, j
+ integer(MPAS_i4) :: matrixSize
+
+ real(MPAS_r8), dimension(:,:), pointer :: matrix
+ real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
+ integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+3 ! extra space for constant vector
+
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,3))
+ allocate(coeffs(matrixSize,3))
+ allocate(pivotIndices(matrixSize))
+
+ matrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
+ coeffs = 0.0_MPAS_r8
+
+ call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+ alpha, matrix, rhs)
+
+ do i = 1, pointCount
+ if(.not. isTangentToInterface(i)) then
+ matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
+ end if
+ matrix(pointCount+1:pointCount+3,i) &
+ = matrix(i,pointCount+1:pointCount+3)
end do
+ do i = 1, 3
+ rhs(pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+ end do
! solve each linear system
- do i = 1, dimensions
+ do i = 1, 3
call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
end do
coefficients = coeffs(1:pointCount,:)
- end subroutine ibInterp_computeVectorNoSlipCoefficientsRBF
+ deallocate(matrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
- !!!!!! Not implemented yet, identical to No Slip !!!!!!!!!
- subroutine ibInterp_computeVectorFreeSlipCoefficientsRBF(pointCount, &
- dimensions, sourcePoints, unitVectors, destinationPoint, &
- alpha, coefficients, dminfo)
+ end subroutine ibInterp_computeVectorFreeSlipRBFCoefficients
- integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
- real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
- real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
- real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
+ subroutine ibInterp_computeVectorFreeSlipPlanarRBFCoefficients(pointCount, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+ alpha, planeBasisVectors, coefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isTangentToInterface
+ integer(MPAS_i4), dimension(pointCount), intent(in) :: normalVectorIndex
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
real(MPAS_r8), intent(in) :: alpha
- real(MPAS_r8), dimension(pointCount, dimensions), intent(out) :: coefficients
+ real(MPAS_r8), dimension(2,3) :: planeBasisVectors
+ real(MPAS_r8), dimension(pointCount, 3), intent(out) :: coefficients
type (dm_info), intent(in) :: dminfo
- integer(MPAS_i4) :: matrixSize, i, j
+ integer(MPAS_i4) :: i, j
+ integer(MPAS_i4) :: matrixSize
+ real(MPAS_r8), dimension(pointCount,2) :: planarSourcePoints
+ real(MPAS_r8), dimension(pointCount,2) :: planarUnitVectors
+ real(MPAS_r8), dimension(2) :: planarDestinationPoint
+
real(MPAS_r8), dimension(:,:), pointer :: matrix
real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
integer(MPAS_i4), dimension(:), pointer :: pivotIndices
- real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+ matrixSize = pointCount+2 ! space for constant vector in plane
- matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,2))
+ allocate(coeffs(matrixSize,2))
+ allocate(pivotIndices(matrixSize))
- allocate(matrix(matrixSize,matrixSize))
- allocate(rhs(matrixSize,dimensions))
- allocate(coeffs(matrixSize,dimensions))
- allocate(pivotIndices(matrixSize))
-
matrix = 0.0_MPAS_r8
rhs = 0.0_MPAS_r8
coeffs = 0.0_MPAS_r8
- do j = 1, pointCount
- do i = 1, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
- rbfValue = evaluateRBF(rSquared)
- matrix(j,i) = rbfValue*unitVectorDotProduct
- end do
- rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
- rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+ do i = 1, pointCount
+ planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:))
+ planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:))
+ planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+ planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
end do
+ planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:))
+ planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
+ call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 2, &
+ planarSourcePoints, isTangentToInterface, normalVectorIndex, planarUnitVectors, &
+ planarDestinationPoint, alpha, matrix, rhs)
do i = 1, pointCount
- matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
- matrix(pointCount+1:pointCount+dimensions,i) &
- = matrix(i,pointCount+1:pointCount+dimensions)
+ if(.not. isTangentToInterface(i)) then
+ matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+ matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+ end if
+ matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
end do
- do i = 1, dimensions
- rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
- end do
+ rhs(pointCount+1,:) = planeBasisVectors(1,:)
+ rhs(pointCount+2,:) = planeBasisVectors(2,:)
! solve each linear system
- do i = 1, dimensions
- call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
- end do
+ call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
+ call LEGS(matrix, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
+
+ coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,1)*coeffs(1:pointCount,2)
+ coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,2)*coeffs(1:pointCount,2)
+ coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,3)*coeffs(1:pointCount,2)
coefficients = coeffs(1:pointCount,:)
- end subroutine ibInterp_computeVectorFreeSlipCoefficientsRBF
+ deallocate(matrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+ end subroutine ibInterp_computeVectorFreeSlipPlanarRBFCoefficients
-
!!!!!!!!!!!!!!!!!!!!!
! private subroutines
!!!!!!!!!!!!!!!!!!!!!
subroutine exit_MPAS(string, dminfo)
- char(len=*), intent(in) :: string
+ character(len=*), intent(in) :: string
type (dm_info), intent(in) :: dminfo
write(0,*) string
@@ -654,96 +758,144 @@
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletMatrix, neumannMatrix, rhs)
- integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+ integer(MPAS_i4), intent(in) :: pointCount
real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
logical, dimension(pointCount), intent(in) :: isInterface
real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
real(MPAS_r8), intent(in) :: alpha
real(MPAS_r8), dimension(pointCount,pointCount), intent(out) :: &
- dirichletCoefficients, neumannCoefficients
+ dirichletMatrix, neumannMatrix
real(MPAS_r8), dimension(pointCount), intent(out) :: rhs
integer(MPAS_i4) :: i, j
-
- integer(MPAS_i4), dimension(matrixSize) :: pivotIndices
-
real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
- dirichletMatrix = 0.0_MPAS_r8
- neumannMatrix = 0.0_MPAS_r8
- rhs = 0.0_MPAS_r8
+ do j = 1, pointCount
+ if(isInterface(j)) then
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ normalDotX = sum(interfaceNormals(j,:) &
+ * (sourcePoints(j,:)-sourcePoints(i,:)))
+ call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+ rbfDerivOverR = rbfDerivOverR/alpha**2
+ dirichletMatrix(i,j) = rbfValue
+ neumannMatrix(i,j) = rbfDerivOverR*normalDotX
+ end do
+ else
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ rbfValue = evaluateRBF(rSquared)
+ dirichletMatrix(i,j) = rbfValue
+ neumannMatrix(i,j) = rbfValue
+ end do
+ end if
+ end do
do j = 1, pointCount
- if(isInterface(j)) cycle
- do i = 1, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- rbfValue = evaluateRBF(rSquared)
- dirichletMatrix(i,j) = rbfValue
- neumannMatrix(i,j) = rbfValue
- end do
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j) = evaluateRBF(rSquared)
end do
+ end subroutine setUpScalarRBFMatrixAndRHS
+
+ subroutine setUpVectorNoSlipRBFMatrixAndRHS(pointCount, dimensions, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, matrix, rhs)
+
+ integer(MPAS_i4), intent(in) :: pointCount, dimensions
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(pointCount,pointCount), intent(out) :: matrix
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(out) :: rhs
+
+ integer(MPAS_i4) :: i, j
+
+ real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+
do j = 1, pointCount
- if(.not. isInterface(j)) cycle
- do i = 1, pointCount
+ do i = j, pointCount
rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- normalDotX = sum(interfaceNormals(j,:) &
- * (sourcePoints(j,:)-sourcePoints(i,:)))
- call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
- rbfDerivOverR = rbfDerivOverR/alpha**2
- dirichletMatrix(i,j) = rbfValue
- neumannMatrix(i,j) = rbfDerivOverR*normalDotX
+ rbfValue = evaluateRBF(rSquared)
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ matrix(i,j) = rbfValue*unitVectorDotProduct
+ matrix(j,i) = matrix(i,j)
end do
end do
do j = 1, pointCount
rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
- rhs(j) = evaluateRBF(rSquared)
+ rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
end do
- end subroutine setUpScalarRBFMatrixAndRHS
+ end subroutine setUpVectorNoSlipRBFMatrixAndRHS
- 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 setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, dimensions, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+ alpha, matrix, rhs)
-!======================================================================
-! BEGINNING OF CROSS_PRODUCT_IN_R3
-!======================================================================
- subroutine cross_product_in_R3(p_1,p_2,p_out)
+ integer(MPAS_i4), intent(in) :: pointCount, dimensions
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isTangentToInterface
+ integer(MPAS_i4), dimension(pointCount), intent(in) :: normalVectorIndex
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(pointCount,pointCount), intent(out) :: matrix
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(out) :: rhs
-!-----------------------------------------------------------------------
-! PURPOSE: compute p_1 cross p_2 and place in p_out
-!-----------------------------------------------------------------------
+ integer(MPAS_i4) :: i, j
-!-----------------------------------------------------------------------
-! intent(in)
-!-----------------------------------------------------------------------
- real (kind=RKIND), intent(in) :: &
- p_1 (3), &
- p_2 (3)
+ real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalVector(dimensions), &
+ normalDotX, unitVectorDotProduct
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
- real (kind=RKIND), intent(out) :: &
- p_out (3)
+ do j = 1, pointCount
+ if(isTangentToInterface(j)) then
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ normalVector = unitVectors(normalVectorIndex(j),:)
+ normalDotX = sum(normalVector * (sourcePoints(j,:)-sourcePoints(i,:)))
+ call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ matrix(i,j) = rbfDerivOverR*unitVectorDotProduct
+ end do
+ else
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ rbfValue = evaluateRBF(rSquared)
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ matrix(i,j) = rbfValue*unitVectorDotProduct
+ end do
+ end if
+ end do
- 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)
+ do j = 1, pointCount
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+ end do
- end subroutine cross_product_in_R3
-!======================================================================
-! END OF CROSS_PRODUCT_IN_R3
-!======================================================================
+ end subroutine setUpVectorFreeSlipRBFMatrixAndRHS
+ 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 cross_product_in_R3(p_1,p_2,p_out)
+ real (kind=RKIND), intent(in) :: p_1 (3), p_2 (3)
+ 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
+
! Updated 10/24/2001.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -948,5 +1100,5 @@
!
END SUBROUTINE ELGS
-end module IB_interpoation
+end module IB_interpolation
</font>
</pre>