diff --git a/.github/scripts/_typos.toml b/.github/scripts/_typos.toml new file mode 100644 index 0000000..b247064 --- /dev/null +++ b/.github/scripts/_typos.toml @@ -0,0 +1,19 @@ +[default] +extend-ignore-identifiers-re = [ + "AttributeID.*Supress.*", +] + +[default.extend-identifiers] +AttributeIDSupressMenu = "AttributeIDSupressMenu" + +[default.extend-words] +iy = "iy" +ths = "ths" +asend = "asend" +alph = "alph" +fo = "fo" +thi = "thi" +doub = "doub" +Doub = "Doub" +Numer = "Numer" +thr = "thr" \ No newline at end of file diff --git a/.github/workflows/pretty.yml b/.github/workflows/pretty.yml index fdfd31e..1931e66 100644 --- a/.github/workflows/pretty.yml +++ b/.github/workflows/pretty.yml @@ -2,9 +2,7 @@ name: Pretty on: push: - pull_request: - workflow_dispatch: jobs: diff --git a/.github/workflows/spelling.yml b/.github/workflows/spelling.yml new file mode 100644 index 0000000..003980d --- /dev/null +++ b/.github/workflows/spelling.yml @@ -0,0 +1,18 @@ +name: Spell Check +on: + push: + pull_request: + workflow_dispatch: + +jobs: + run: + name: Spell Check + runs-on: ubuntu-latest + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Spell Check + uses: crate-ci/typos@master + with: + config: ${{github.workspace}}/.github/scripts/_typos.toml \ No newline at end of file diff --git a/.gitignore b/.gitignore index 4f18b65..4cbc252 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ */*/v.* */*/D/* !*/*/D/.gittouch +*/*/*/*.log *.o *.mod .depend @@ -11,4 +12,5 @@ libcommon.a .vscode field packages -traction \ No newline at end of file +traction + diff --git a/Makefile.in b/Makefile.in index ca00b91..0eb9be4 100755 --- a/Makefile.in +++ b/Makefile.in @@ -1,10 +1,12 @@ -# Directories (change this to where you put RBC3D!) -WORK_DIR = /storage/home/hcoda1/6/sbryngelson3/p-sbryngelson3-0/test-rbc/RBC3D +# Get RBC3D root directory +WORK_DIR := $(shell dirname $(realpath $(lastword $(MAKEFILE_LIST)))) # Package directories (shouldn't need to change these) -PETSC_DIR = $(WORK_DIR)/packages/mypetsc -include $(PETSC_DIR)/conf/variables +PETSC_DIR = $(WORK_DIR)/packages/petsc-3.19.6 +PETSC_ARCH = $(PETSC_DIR)/petsc_configure + SPHEREPACK_DIR = $(WORK_DIR)/packages/spherepack3.2 + LAPACK_DIR = $(WORK_DIR)/packages/lapack-3.11 # Makedependf90 binary @@ -18,19 +20,21 @@ vpath $(WORK_DIR)/common # Includes PETSC_INCLUDE = $(PETSC_DIR)/include +PETSC_ARCH_INCLUDE = $(PETSC_ARCH)/include NETCDF_INCLUDE = $(NETCDF_DIR)/include COMMON_INCLUDE = -I$(WORK_DIR)/common -INCLUDE = $(COMMON_INCLUDE) -I$(PETSC_INCLUDE) -I$(NETCDF_INCLUDE) +INCLUDE = $(COMMON_INCLUDE) -I$(PETSC_INCLUDE) -I$(PETSC_ARCH_INCLUDE) -I$(NETCDF_INCLUDE) # Libraries COMMON_LIB = $(WORK_DIR)/common/libcommon.a SPHPK_LIB = -L$(SPHEREPACK_DIR)/lib -lspherepack FFTW_LIB = -L$(FFTW_DIR)/lib -lfftw3 NETCDF_LIB = -L$(NETCDF_DIR)/lib -lnetcdff -PETSC_LIB = -L$(PETSC_DIR)/lib $(PETSC_KSP_LIB_BASIC) -MKL_LIB = -L$(MKL_ROOT)lib/intel64/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl -LAPACK_LIB = -L$(LAPACK_DIR) -llapack +PETSC_LIB = -Wl,-rpath,$(PETSC_ARCH)/lib -L$(PETSC_ARCH)/lib -lpetsc -lstdc++ +MKL_LIB = -L$(MKL_ROOT)lib/intel64/ -lmkl_gf_lp64 -lmkl_core -lmkl_sequential -lpthread -lm -ldl +BLAS_LIB = -L$(BLAS_DIR) -lrefblas +LAPACK_LIB = -L$(LAPACK_DIR) -llapack -lrefblas # Compiler and linker FC = mpif90 diff --git a/README.md b/README.md index 9a7fee1..6abe811 100644 --- a/README.md +++ b/README.md @@ -23,8 +23,22 @@ This codebase solves the boundary integral form of the Stokes equations via an a ### Installation -Building RBC3D can be a fragile process and we are working on improving it. -In the meantime, a careful documentation of build instructions [is available here](https://github.com/comp-physics/RBC3D/blob/master/install/readme.md). +To install on PACE Phoenix, you can salloc a node and then run: + +```shell +bash rbc.sh install-with-mkl +``` + +Then to execute and run a case, you can: +```shell +cd examples/case +make .depend +make +srun -n 1 ./initcond +srun ./tube +``` + +On other supercomputing clusters, it should be easy to replace line of 7 `./install/install-with-mkl.sh` with the modules available on your system and change the directories in `Makefile.in` to point to those modules. If one of these isn't available, you can follow the manual build instructions [available here](https://github.com/comp-physics/RBC3D/blob/master/install/readme.md). ### Papers that use RBC3D diff --git a/common/ModBasicMath.F90 b/common/ModBasicMath.F90 index 208c107..d7f590d 100644 --- a/common/ModBasicMath.F90 +++ b/common/ModBasicMath.F90 @@ -1,4 +1,4 @@ -! Collection of basic math operation +! Collection of basic math operations module ModBasicMath use ModDataTypes @@ -7,7 +7,7 @@ module ModBasicMath private - public :: VecNorm, & + public :: VecNormL2, & CrossProd, & InvMat2, & InvMat3, & @@ -27,12 +27,12 @@ module ModBasicMath !********************************************************************** ! L2 norm of a vector ! |a| - function VecNorm(a) result(c) + function VecNormL2(a) result(c) real(WP) :: a(:), c c = sqrt(sum(a*a)) - end function VecNorm + end function VecNormL2 !********************************************************************** ! Cross product of two vectors @@ -167,7 +167,9 @@ subroutine Matrix_PseudoInvert(A, B) AtA = matmul(B, A) ! Solve (AtA)X = B and store solution in B - call DPOSV('U', N, M, Ata, N, B, N, ierr) + call dposv('U', N, M, Ata, N, B, N, ierr) + + if (ierr .ne. 0) write (*, *) "Matrix_PseudoInvert error code: ", ierr ! Deallocate working arrays deallocate (AtA) @@ -212,7 +214,9 @@ subroutine QuadFit_1D(x, f, a0, a1, a2) lhs(3, 1) = lhs(1, 3) lhs(3, 2) = lhs(2, 3) - call DPOSV('U', 3, 1, lhs, 3, rhs, 3, ierr) + call dposv('U', 3, 1, lhs, 3, rhs, 3, ierr) + + if (ierr .ne. 0) write (*, *) "QuadFit_1D error code: ", ierr a0 = rhs(1) a1 = rhs(2) @@ -222,6 +226,7 @@ end subroutine QuadFit_1D !********************************************************************** ! 2D Quadratic fit +! Uses least squares normal equation method ! Arguments: ! x(i,:), f(i) -- coordinates and functional value of the i-th point ! Note: @@ -265,7 +270,9 @@ subroutine QuadFit_2D(x, f, a0, a1, a2, a11, a12, a22) end do ! jj end do ! ii - call DPOSV('U', 6, 1, lhs, 6, rhs, 6, ierr) + call dposv('U', 6, 1, lhs, 6, rhs, 6, ierr) + + if (ierr .ne. 0) write (*, *) "QuadFit_2D error code: ", ierr a0 = rhs(1) a1 = rhs(2) @@ -410,7 +417,7 @@ subroutine BsplineFunc(xc, P, imin, w) end subroutine BsplineFunc !********************************************************************** -! Random number generator, copied from Numerical Recipies +! Random number generator, copied from Numerical Recipes function RandomNumber(idum) integer, intent(IN), optional :: idum real(WP) :: RandomNumber diff --git a/common/ModConf.F90 b/common/ModConf.F90 index 2ee3b4e..439590f 100644 --- a/common/ModConf.F90 +++ b/common/ModConf.F90 @@ -5,6 +5,9 @@ module ModConf use ModDataTypes use ModDataStruct +#include "petsc/finclude/petsc.h" + use petsc + implicit none ! Domain size @@ -108,7 +111,6 @@ module ModConf subroutine InitMPI(split_comm) logical, optional :: split_comm -#include "../petsc_include.h" integer :: numNodes, nodeNum character(MPI_Max_Processor_Name) :: machinename integer :: lenname @@ -207,8 +209,6 @@ end subroutine InitMPI !********************************************************************** ! Finalize MPI and PETSc subroutine FinalizeMPI - -#include "../petsc_include.h" integer :: ierr if (MPI_COMM_Ewald /= MPI_COMM_WORLD) then @@ -437,7 +437,7 @@ subroutine DomainDecomp end subroutine DomainDecomp !********************************************************************** -! Wheather a source point make contribution to the local domain +! Whether a source point make contribution to the local domain function Is_Source(x) real(WP) :: x(3) logical :: Is_Source @@ -493,7 +493,7 @@ function Cell_Has_Source(cell) result(hasSource) end function Cell_Has_Source !********************************************************************** -! Wheather a triangle has active source points +! Whether a triangle has active source points ! Arguments: ! x(i,:) -- the coordinates of the i-th vertex function Tri_Has_Source(x) result(hasSource) diff --git a/common/ModDataStruct.F90 b/common/ModDataStruct.F90 index 3b36a91..90def19 100644 --- a/common/ModDataStruct.F90 +++ b/common/ModDataStruct.F90 @@ -2,8 +2,10 @@ module ModDataStruct use ModDataTypes +#include "petsc/finclude/petsc.h" + use petsc + implicit none -#include "../petsc_include.h" private @@ -140,8 +142,8 @@ module ModDataStruct ! radius -- patch radius ! nrad, nazm -- number of patch points along radial and ! azimuthal directions -! thL, phiL -- local coordiantes of patch mesh points -! w -- weight of integration weight (including masking fucntion) +! thL, phiL -- local coordinates of patch mesh points +! w -- weight of integration weight (including masking function) ! thG(:,:,i,j), phiG(:,:,i,j) -- global latitudinal and longitudinal ! coordinates of a patch centered at (i,j)-th mesh point type t_RbcPolarPatch @@ -172,7 +174,7 @@ module ModDataStruct ! ! a3 -- element normal ! area -- element area -! epsDist -- threshhold distance for singular integration +! epsDist -- threshold distance for singular integration ! areaTot -- total surface area ! ! lhs -- lhs matrix for no-slip boundary condition @@ -198,7 +200,7 @@ module ModDataStruct ! nPoint -- number of points ! x(i,:) -- coordinate of the i-th point ! f(i,:) -- force singularity at the i-th point -! g(i,:) -- double-layer potential source singularity at the i-th piont +! g(i,:) -- double-layer potential source singularity at the i-th point ! a3(i,:) -- normal direction ! lam(i) -- viscRatio ! @@ -207,7 +209,7 @@ module ModDataStruct ! for wall surface, (surfId, iele, -1) ! ! Nc -- number of cells -! iLbNc -- Nc/Lb, for indentifying which cell a point lies in +! iLbNc -- Nc/Lb, for identifying which cell a point lies in ! hoc -- first point in the cell ! next -- next point type t_SourceList @@ -227,7 +229,7 @@ module ModDataStruct ! nPoint -- number of points ! x -- coordinate and velocities of points ! indx -- same as those in t_SourceList -! active -- wheather the target point is active +! active -- Whether the target point is active type t_TargetList integer :: nPoint real(WP), pointer, dimension(:, :) :: x diff --git a/common/ModHashTable.F90 b/common/ModHashTable.F90 index 8e4c021..33bd9fc 100644 --- a/common/ModHashTable.F90 +++ b/common/ModHashTable.F90 @@ -60,7 +60,7 @@ end subroutine HashTable_Build !********************************************************************** ! Compute the Hash table index of a point ! Argument: -! Nc -- nunber of cell blocks +! Nc -- number of cell blocks ! iLbNc -- Nc/(local domain size) ! x -- point coordinate ! i1, i2, i3 -- Hash index of x(:) diff --git a/common/ModIO.F90 b/common/ModIO.F90 index ad72e6d..47aa121 100644 --- a/common/ModIO.F90 +++ b/common/ModIO.F90 @@ -210,9 +210,9 @@ subroutine WriteManyRBCs(fn, nrbc, rbcs) write (cell_unit, '(A,I9,A,I9,A)') 'ZONE I=', nlat + 1, ' J=', nlon + 1, ' F=POINT' do ilon = 1, nlon - do ilat = 0, nlat - write (cell_unit, '(3F20.10)') x(ilat, ilon, :) - end do ! ilat + do ilat = 0, nlat + write (cell_unit, '(3F20.10)') x(ilat, ilon, :) + end do ! ilat end do ! ilon do ilat = 0, nlat write (cell_unit, '(3F20.10)') x(ilat, 1, :) @@ -230,7 +230,7 @@ end subroutine WriteManyRBCs ! Write the shape of blood cells of specified type to file ! Arguments: ! fn -- file suffix name -! nrbc -- nubmer of cells +! nrbc -- number of cells ! rbcs -- blood cells ! type -- type filter (1: rbc, 2: leukocyte, 3: sickle cell) subroutine WriteManyRBCsByType(fn, nrbc, rbcs, type) diff --git a/common/ModIntOnWalls.F90 b/common/ModIntOnWalls.F90 index 54ac784..6d417b7 100644 --- a/common/ModIntOnWalls.F90 +++ b/common/ModIntOnWalls.F90 @@ -10,6 +10,9 @@ module ModIntOnWalls use ModBasicMath use ModEwaldFunc +#include "petsc/finclude/petsc.h" + use petsc + implicit none private @@ -61,7 +64,7 @@ subroutine AddIntOnWalls(c1, tlist, v) call SingIntOnWall(c1, wall, vtmp) do i = 1, wall%nvert v(p + i, :) = v(p + i, :) + vtmp(i, :)/tlist%Acoef(p + i) !COEF -! v(p+i,:) = v(p+i,:) + vtmp(i,:)/(1.+tlist%lam(p+i)) !COEF + ! v(p+i,:) = v(p+i,:) + vtmp(i,:)/(1.+tlist%lam(p+i)) !COEF end do deallocate (vtmp) @@ -134,34 +137,43 @@ subroutine SingIntOnWall(c1, wall, v) real(WP) :: c1 type(t_Wall) :: wall real(WP) :: v(:, :) + real(WP), allocatable :: v1D(:), f1D(:) -#include "../petsc_include.h" - integer :: nrow, i + integer :: nrow, i, j integer, allocatable :: irows(:) Vec :: f_vec, lhsf_vec integer :: ierr + integer :: Mat_m, Mat_n nrow = 3*wall%nvert + allocate (irows(nrow)) + allocate (f1D(nrow)) + allocate (v1D(nrow)) + + f1D = reshape(wall%f, (/nrow/)) + call VecCreateSeq(PETSC_COMM_SELF, nrow, f_vec, ierr) call VecCreateSeq(PETSC_COMM_SELF, nrow, lhsf_vec, ierr) irows = (/(i, i=0, nrow - 1)/) - call VecSetValues(f_vec, nrow, irows, wall%f, INSERT_VALUES, ierr) + call VecSetValues(f_vec, nrow, irows, f1D, INSERT_VALUES, ierr) call VecAssemblyBegin(f_vec, ierr) + call MatMult(wall%lhs, f_vec, lhsf_vec, ierr) - call VecGetValues(lhsf_vec, nrow, irows, v, ierr) + call VecGetValues(lhsf_vec, nrow, irows, v1D, ierr) + + v = reshape(v1D, (/wall%nvert, 3/)) v = c1*v call VecDestroy(f_vec, ierr) call vecDestroy(lhsf_vec, ierr) - deallocate (irows) - + deallocate (irows, f1D, v1D) end subroutine SingIntOnWall !********************************************************************** -! Prepare for the self-intergration on the wall +! Prepare for the self-integration on the wall ! Arguments: ! wall -- ! @@ -170,7 +182,6 @@ end subroutine SingIntOnWall subroutine PrepareSingIntOnWall(wall) type(t_wall) :: wall -#include "../petsc_include.h" type(t_targetlist), pointer :: tlist type(t_sourcelist), pointer :: slist integer :: nvert, nrow @@ -232,7 +243,8 @@ subroutine PrepareSingIntOnWall(wall) end do ! j1 end do ! i - nnz = nnz*2 ! a very conservative estimate + ! nnz = nnz*2 ! a very conservative estimate + nnz = nnz*12 ! a very conservative estimate call MatCreateSeqAIJ(PETSC_COMM_SELF, nrow, nrow, 0, nnz, wall%lhs, ierr) ! Assemble the matrix @@ -279,7 +291,7 @@ subroutine PrepareSingIntOnWall(wall) ivert = wall%e2v(iele, l) icols = (/ivert, ivert + nvert, ivert + 2*nvert/) - 1 values = reshape(transpose(lhs(l, :, :)), (/9/)) - call MatSetValues(wall%lhs, 3, irows, 3, icols, values, ADD_VAlUES, ierr) + call MatSetValues(wall%lhs, 3, irows, 3, icols, values, ADD_VALUES, ierr) end do ! l 999 j = slist%next(j) @@ -299,7 +311,7 @@ end subroutine PrepareSingIntOnWall !********************************************************************** ! Compute the regular surface integral over a triangle ! Arguments: -! x(i,:), f(i,:) -- trianglular vertex coordinates and force densities +! x(i,:), f(i,:) -- triangular vertex coordinates and force densities ! xtar -- target point ! rhs -- integration ! lhs -- influence matrix @@ -354,7 +366,7 @@ end subroutine Tri_Int_Regular !********************************************************************** ! Use Duffy's rule to compute the surface integral over a triangle ! Arguments: -! x(i,:), f(i,:) -- trianglular vertex coordinates and force densities +! x(i,:), f(i,:) -- triangular vertex coordinates and force densities ! xtar -- target point ! (s0, t0) -- reference coordinate of the nearest point to xtar ! rhs -- rhs diff --git a/common/ModNoSlip.F90 b/common/ModNoSlip.F90 index 8f2c581..88e9cdd 100644 --- a/common/ModNoSlip.F90 +++ b/common/ModNoSlip.F90 @@ -13,9 +13,10 @@ module ModNoSlip use ModTargetList use ModSphPk - implicit none +#include "petsc/finclude/petsc.h" + use petsc -#include "../petsc_include.h" + implicit none ! npoint -- total number of wall mesh points, with redundancy ! dof -- degrees of freedom @@ -56,7 +57,7 @@ subroutine NoSlipWall ! Initialize the solver if (.not. solver_inited) then - ! Calcuate system size + ! Calculate system size npoint = 0 dof = 0 do iwall = 1, nwall @@ -73,7 +74,8 @@ subroutine NoSlipWall call MatShellSetOperation(mat_lhs, MATOP_MULT, MyMatMult, ierr) call KSPCreate(PETSC_COMM_SELF, ksp_lhs, ierr) - call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, SAME_NONZERO_PATTERN, ierr) + ! call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, SAME_NONZERO_PATTERN, ierr) + call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, ierr) call KSPSetType(ksp_lhs, KSPGMRES, ierr) call KSPGetPC(ksp_lhs, pc_lhs, ierr) ! call KSPSetInitialGuessNonzero(ksp_lhs, PETSC_TRUE, ierr) ! TRIAL @@ -81,7 +83,7 @@ subroutine NoSlipWall !SHB modify from 1.D-3 to 1.D-6 and now to eps_Ewd call KSPSetTolerances(ksp_lhs, eps_Ewd, & - PETSC_DEFAULT_DOUBLE_PRECISION, PETSC_DEFAULT_DOUBLE_PRECISION, & + PETSC_DEFAULT_REAL, PETSC_DEFAULT_REAL, & 60, ierr) !!$ print * !!$ print *, "SMALL WALL ITERATIONS" diff --git a/common/ModPFFTW.F90 b/common/ModPFFTW.F90 index d8c64dd..1bca14d 100644 --- a/common/ModPFFTW.F90 +++ b/common/ModPFFTW.F90 @@ -68,7 +68,7 @@ SUBROUTINE Init_PFFTW(ixBgn, ixEnd, iqBgn, iqEnd) end do end if - pointsR(0, 1) = 1 !indicies of R on processor + pointsR(0, 1) = 1 !indices of R on processor pointsR(0, 2) = chunkR(0) pointsT(0, 1) = 1 pointsT(0, 2) = chunkT(0) diff --git a/common/ModPME.F90 b/common/ModPME.F90 index 0b77e25..70324c0 100644 --- a/common/ModPME.F90 +++ b/common/ModPME.F90 @@ -387,7 +387,7 @@ subroutine Update_Buff_Vel call MPI_SendRecv(bufSend, sizebuf, MPI_WP, nodeRight, nodeNum, & bufRecv, sizeBuf, MPI_WP, nodeLeft, nodeLeft, & MPI_COMM_Ewald, stat, ierr) -! print*,' callled mpi' +! print*,' called mpi' vv(:, :, ixBgn(3) - PBspln:ixBgn(3) - 1, :) = reshape(bufRecv, (/Nb(1), Nb(2), PBspln, 3/)) ! Deallocate working arrays diff --git a/common/ModPolarPatch.F90 b/common/ModPolarPatch.F90 index 6420490..add12a2 100644 --- a/common/ModPolarPatch.F90 +++ b/common/ModPolarPatch.F90 @@ -148,7 +148,7 @@ subroutine PolarPatch_Build(th0, phi0, thL, phiL, thG, phiG) end subroutine PolarPatch_Build !********************************************************************** -! Find the Cartesian mesh pionts within a polar patch +! Find the Cartesian mesh points within a polar patch ! Arguments: ! (th0, phi0) -- origin of the patch ! r0 -- radius of the patch diff --git a/common/ModPolyRoots.F90 b/common/ModPolyRoots.F90 index 78d494e..ccab257 100644 --- a/common/ModPolyRoots.F90 +++ b/common/ModPolyRoots.F90 @@ -2,7 +2,7 @@ MODULE ModPolyRoots ! --------------------------------------------------------------------------- ! PURPOSE - Solve for the roots of a polynomial equation with real -! coefficients, up to quartic order. Retrun a code indicating the nature +! coefficients, up to quartic order. Return a code indicating the nature ! of the roots found. ! AUTHORS - Alfred H. Morris, Naval Surface Weapons Center, Dahlgren,VA @@ -387,7 +387,7 @@ SUBROUTINE QuarticRoots(a, z) ! SOLVE THE RESOLVENT CUBIC EQUATION. THE CUBIC HAS AT LEAST ONE ! NONNEGATIVE REAL ROOT. IF W1, W2, W3 ARE THE ROOTS OF THE CUBIC -! THEN THE ROOTS OF THE ORIGINIAL EQUATION ARE +! THEN THE ROOTS OF THE ORIGINAL EQUATION ARE ! Z = -B + CSQRT(W1) + CSQRT(W2) + CSQRT(W3) ! WHERE THE SIGNS OF THE SQUARE ROOTS ARE CHOSEN SO ! THAT CSQRT(W1) * CSQRT(W2) * CSQRT(W3) = -Q/8. diff --git a/common/ModPostProcess.F90 b/common/ModPostProcess.F90 index c0b029c..9e091e9 100644 --- a/common/ModPostProcess.F90 +++ b/common/ModPostProcess.F90 @@ -22,7 +22,7 @@ module ModPostProcess contains !********************************************************************** -! Compute the velocity on an abitrary point sets +! Compute the velocity on an arbitrary point sets ! Arguments: ! tlist -- the target list ! v(:,:) -- velocities @@ -167,7 +167,7 @@ function DistFromWall(type) result(minDist) do ivert = 1, wall%nvert do ilat = 1, rbc%nlat do ilon = 1, rbc%nlon - currDist = VecNorm((rbc%x(ilat, ilon, :)) - (wall%x(ivert, :))) + currDist = NORM2((rbc%x(ilat, ilon, :)) - (wall%x(ivert, :))) if (currDist .le. minDist) then minDist = currDist cellPoint = rbc%x(ilat, ilon, :) diff --git a/common/ModQuadRule.F90 b/common/ModQuadRule.F90 index 6396821..545d0b8 100644 --- a/common/ModQuadRule.F90 +++ b/common/ModQuadRule.F90 @@ -131,7 +131,7 @@ end subroutine GaussQuad_Finalize ! double precision. ! ! Note: -! Copied from "Numerical Recipies in Fortran 90" +! Copied from "Numerical Recipes in Fortran 90" subroutine GauLeg(x1, x2, n, x, w) real(WP) :: x1, x2 integer :: n @@ -212,7 +212,7 @@ end subroutine GauLeg ! Arguments: ! xmin, xmax -- integration interval ! (a, b) -- the target point coordinate. b is very small (but non-zero), -! which will cause a near singularity for the integration arond +! which will cause a near singularity for the integration around ! (a,0) ! n, x, w -- same as those in subroutine GauLeg ! diff --git a/common/ModRbc.F90 b/common/ModRbc.F90 index 161be95..5a2bb41 100644 --- a/common/ModRbc.F90 +++ b/common/ModRbc.F90 @@ -453,7 +453,7 @@ subroutine RBC_ComputeGeometry(cell) cell%detj(ilat, ilon) = sqrt(detA)/sin(cell%th(ilat)) a3 = CrossProd(a1, a2) - cell%a3(ilat, ilon, :) = a3/VecNorm(a3) + cell%a3(ilat, ilon, :) = a3/NORM2(a3) end do ! ilat end do ! ilon @@ -509,7 +509,7 @@ subroutine RBC_ComputeGeometry(cell) end subroutine RBC_ComputeGeometry !********************************************************************** -! Compute the covariant gradient of a vecotr field +! Compute the covariant gradient of a vector field ! Arguments: ! cell -- ! v(ilat,ilon,:) -- covariant vector component @@ -921,7 +921,7 @@ subroutine Shell_ElasStrs(cell, cellRef, tau) do ilat = 1, nlat do ilon = 1, nlon - ! V2 is the strech tensor + ! V2 is the stretch tensor V2 = matmul(cell%a(ilat, ilon, :, :), cellRef%a_rcp(ilat, ilon, :, :)) lbd1 = V2(1, 1) + V2(2, 2) - 2.0 diff --git a/common/ModRbcSingInt.F90 b/common/ModRbcSingInt.F90 index 847c1c5..2760324 100644 --- a/common/ModRbcSingInt.F90 +++ b/common/ModRbcSingInt.F90 @@ -23,7 +23,7 @@ module ModRbcSingInt contains !********************************************************************** -! Comptue the correction due to singular integration +! Compute the correction due to singular integration ! Arguments: ! dv -- correction to the simple point-point sum subroutine RBC_SingInt(c1, c2, rbc, ilat0, ilon0, dv) @@ -45,9 +45,9 @@ subroutine RBC_SingInt(c1, c2, rbc, ilat0, ilon0, dv) ! Add integrals on the local patch ! Note: ! The spline interpolation of surface coordinates does not exactly - ! coincide with the orignal ones on mesh points. + ! coincide with the original ones on mesh points. ! - ! Since the coordinates of quadrature poins are by interpolation, we + ! Since the coordinates of quadrature points are by interpolation, we ! should also use interpolation for the target point position, even ! though the target point. ! @@ -187,7 +187,7 @@ subroutine RBC_NearSingInt_Subtract(C1, C2, rbc, xi, x0, th0, phi0, radPat, dv) ! Initialize dv = 0. - ! Substract the imporperly added contribution from the surface + ! subtract the imporperly added contribution from the surface allocate (ijsPat(rbc%nlat*rbc%nlon, 2)) call PolarPatch_FindPoints(th0, phi0, radPat, rbc%th, rbc%phi, nptPat, ijsPat) diff --git a/common/ModRepulsion.F90 b/common/ModRepulsion.F90 index ee361cb..e882df4 100644 --- a/common/ModRepulsion.F90 +++ b/common/ModRepulsion.F90 @@ -165,7 +165,7 @@ subroutine AddIntraCellRepulsionForce thi = rbc%th(ilat) phii = rbc%phi(ilon) - ! Set distance threshhold + ! Set distance threshold epsDist = 2*sqrt(rbc%area)/rbc%nlat epsDistRef = rbc%patch%radius @@ -234,7 +234,7 @@ subroutine AddIntraCellRepulsionForce end do ! i end if - ! Output informaion + ! Output information call MPI_AllReduce(distMin, distMin_Glb, 1, MPI_WP, MPI_MIN, MPI_COMM_WORLD, ierr) call MPI_AllReduce(cntDf, cntDf_Glb, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr) @@ -304,7 +304,7 @@ subroutine InterCellRepulsion xi = tlist%x(i, :) surfId_i = tlist%indx(i, 0) - ! Set up distance threshhold + ! Set up distance threshold rbc => rbcs(surfId_i) call Closest_Neighbor_Cell(xi, surfId_i, epsDist, x1, dist1) @@ -474,7 +474,7 @@ end subroutine LeukWallRepulsion ! Argument: ! xi -- the point ! surfId_i -- ID of the surface that xi(:) lies on -! epsDist -- distance threshhold +! epsDist -- distance threshold ! x0 -- projection of x0 on the closest neighbor surface ! dist0 -- ||xi - x0|| subroutine Closest_Neighbor_Cell(xi, surfId_i, epsDist, x0, dist0) @@ -550,7 +550,7 @@ end subroutine Closest_Neighbor_Cell ! Arguments: ! xi -- the point ! surfId_i -- ID of the surface that xi lies on -! epsDist -- distance threshhold +! epsDist -- distance threshold ! x0 -- the projection of xi on the closest neighboring wall ! dist0 -- ||xi -x0|| subroutine Closest_Neighbor_Wall(xi, surfId_i, epsDist, x0, dist0) diff --git a/common/ModSpline.F90 b/common/ModSpline.F90 index 6660bc6..230bbef 100644 --- a/common/ModSpline.F90 +++ b/common/ModSpline.F90 @@ -206,7 +206,7 @@ subroutine Spline_FindProjection(spln, xtar, th0, phi0, x0) integer :: iter integer, parameter :: iterMax = 3 - integer, parameter :: nth = 2 ! number of points on the patch arond (th0, phi0) + integer, parameter :: nth = 2 ! number of points on the patch around (th0, phi0) integer, parameter :: nphi = 8 real(WP) :: h real(WP) :: thPat_L(nth), phiPat_L(nphi), thPat(nth, nphi), phiPat(nth, nphi) diff --git a/common/ModTimeInt.F90 b/common/ModTimeInt.F90 index 1da24bc..60e2bb8 100644 --- a/common/ModTimeInt.F90 +++ b/common/ModTimeInt.F90 @@ -398,7 +398,7 @@ subroutine TimeIntModVC time = time0 do lt = Nt0 + 1, Nt - clockBgn = MPI_WTime() ! Start timeing + clockBgn = MPI_WTime() ! Start timing ! Evolve cells ! print *,"NO VEL" @@ -580,7 +580,7 @@ subroutine OneTimeInt time = time0 lt = Nt0 + 1 - clockBgn = MPI_WTime() ! Start timeing + clockBgn = MPI_WTime() ! Start timing ! rbc => rbcs(1) ! allocate(xsave(rbc%nlat,rbc%nlon,3)) @@ -667,7 +667,7 @@ subroutine OneTimeIntModVC time = time0 lt = Nt0 + 1 - clockBgn = MPI_WTime() ! Start timeing + clockBgn = MPI_WTime() ! Start timing ! Evolve cells ! print *,"NO VEL" diff --git a/common/ModVelSolver.F90 b/common/ModVelSolver.F90 index 04023a7..e3dbfaa 100644 --- a/common/ModVelSolver.F90 +++ b/common/ModVelSolver.F90 @@ -13,9 +13,10 @@ module ModVelSolver use ModTargetList use ModSphPk - implicit none +#include "petsc/finclude/petsc.h" + use petsc -#include "../petsc_include.h" + implicit none Mat :: mat_lhs Vec :: vec_rhs, vec_sol @@ -58,7 +59,7 @@ subroutine Solve_RBC_Vel(v) ! Initialize the solver if (.not. solver_inited) then - ! Calcuate matrix size + ! Calculate matrix size dof = 0 do irbc = 1, nrbc rbc => rbcs(irbc) @@ -77,15 +78,16 @@ subroutine Solve_RBC_Vel(v) call MatShellSetOperation(mat_lhs, MATOP_MULT, MyMatMult, ierr) call KSPCreate(PETSC_COMM_SELF, ksp_lhs, ierr) - call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, SAME_NONZERO_PATTERN, ierr) + ! call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, SAME_NONZERO_PATTERN, ierr) + call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, ierr) call KSPSetType(ksp_lhs, KSPGMRES, ierr) call KSPGetPC(ksp_lhs, pc_lhs, ierr) call KSPSetInitialGuessNonzero(ksp_lhs, PETSC_TRUE, ierr) ! TRIAL call PCSetType(pc_lhs, PCNONE, ierr) call KSPSetTolerances(ksp_lhs, 1.D-11, & - PETSC_DEFAULT_DOUBLE_PRECISION, & - PETSC_DEFAULT_DOUBLE_PRECISION, & + PETSC_DEFAULT_REAL, & + PETSC_DEFAULT_REAL, & 200, ierr); print *, "SETTING HIGH MAX ITER/ERR" solver_inited = .true. @@ -382,7 +384,7 @@ subroutine DeflateGetSBModes ! w \cross (x-xc) ! ! All are includeded in the GS orthonomalization for simplicity... - ! ... it would be straighforward to normalize the translation + ! ... it would be straightforward to normalize the translation do ilon = 1, nlon do ilat = 1, nlat @@ -635,7 +637,7 @@ end subroutine Proj ! ! Important: ! For the spherical harmonic coefficiets array, the first dimension -! is for longitudinal direciton, while the second is for the colatitudinal +! is for longitudinal direction, while the second is for the colatitudinal ! direction. subroutine Glob_Sph_Trans(v, c, direction) real(WP) :: v(:, :), c(:) diff --git a/common/ModWall.F90 b/common/ModWall.F90 index 02b1db2..9e76c92 100644 --- a/common/ModWall.F90 +++ b/common/ModWall.F90 @@ -20,7 +20,7 @@ module ModWall ! Create a wall ! Arguments: ! wall -- -! nvert, nele -- number of vertices andd elements +! nvert, nele -- number of vertices and elements subroutine Wall_Create(wall, nvert, nele) type(t_wall) :: wall integer :: nvert, nele diff --git a/examples/Makefile.lapack b/examples/Makefile.lapack new file mode 100644 index 0000000..b5c3cb9 --- /dev/null +++ b/examples/Makefile.lapack @@ -0,0 +1,22 @@ +include ../../Makefile.in + +LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(BLAS_LIB) $(LAPACK_LIB) $(STATIC) +EXECUTABLES = tube initcond + +all : $(EXECUTABLES) + +lib : $(wildcard $(WORK_DIR)/common/*.F90) + make -C $(WORK_DIR)/common + +$(EXECUTABLES) : % : %.o $(wildcard $(COMMON_DIR)/*.F90) + make lib + $(FC) $(LDFLAGS) -o $@ $< $(LIB) + +clean : + -rm -f $(EXECUTABLES) *.o *.mod *.a core + +# Dependency +.depend : $(wildcard *.F90) + $(MAKEDEPEND_BIN) $^ > .depend + +include .depend diff --git a/examples/Makefile.mkl b/examples/Makefile.mkl new file mode 100644 index 0000000..51763b3 --- /dev/null +++ b/examples/Makefile.mkl @@ -0,0 +1,22 @@ +include ../../Makefile.in + +LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC) +EXECUTABLES = tube initcond + +all : $(EXECUTABLES) + +lib : $(wildcard $(WORK_DIR)/common/*.F90) + make -C $(WORK_DIR)/common + +$(EXECUTABLES) : % : %.o $(wildcard $(COMMON_DIR)/*.F90) + make lib + $(FC) $(LDFLAGS) -o $@ $< $(LIB) + +clean : + -rm -f $(EXECUTABLES) *.o *.mod *.a core + +# Dependency +.depend : $(wildcard *.F90) + $(MAKEDEPEND_BIN) $^ > .depend + +include .depend diff --git a/examples/carotid_web/Makefile b/examples/carotid_web/Makefile index 149406a..c68a502 100755 --- a/examples/carotid_web/Makefile +++ b/examples/carotid_web/Makefile @@ -1,22 +1 @@ -include ../../Makefile.in - -LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(LAPACK_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC) -EXECUTABLES = initcond tube - -all : $(EXECUTABLES) - -lib : $(wildcard $(WORK_DIR)/common/*.F90) - make -C $(WORK_DIR)/common - -$(EXECUTABLES) : % : %.o $(wildcard $(COMMON_DIR)/*.F90) - make lib - $(FC) $(LDFLAGS) -o $@ $< $(LIB) - -clean : - -rm -f $(EXECUTABLES) *.o *.mod *.a core - -# Dependency -.depend : $(wildcard *.F90) - $(MAKEDEPEND_BIN) $^ > .depend - -include .depend +include ../Makefile.mkl \ No newline at end of file diff --git a/examples/carotid_web/initcond.F90 b/examples/carotid_web/initcond.F90 index 040dd7f..83d7295 100644 --- a/examples/carotid_web/initcond.F90 +++ b/examples/carotid_web/initcond.F90 @@ -326,7 +326,7 @@ logical function check_wall_collision(cell) sp = walls(i2)%x(j2, :) - if (VecNorm(c1p - sp) .lt. threshold) then + if (NORM2(c1p - sp) .lt. threshold) then check_wall_collision = .true. return end if @@ -353,7 +353,7 @@ logical function check_cell_collision(cell1, cell2) check_cell_collision = .false. !first check the centers; if distance > 4 then we are sure that they don't collide - if (VecNorm(cell1%xc - cell2%xc) .ge. 4) return + if (NORM2(cell1%xc - cell2%xc) .ge. 4) return !each point in cell1 do i = 1, cell1%nlat diff --git a/examples/carotid_web/tube.F90 b/examples/carotid_web/tube.F90 index 9e2f7c3..ccbc49c 100755 --- a/examples/carotid_web/tube.F90 +++ b/examples/carotid_web/tube.F90 @@ -12,12 +12,13 @@ program cells_in_tube use ModBasicMath use ModQuadRule +#include "petsc/finclude/petsc.h" + use petsc + implicit none integer :: cutoff character(CHRLEN) :: fn -#include "../../petsc_include.h" - call InitAll call TimeInt_Euler @@ -31,7 +32,7 @@ program cells_in_tube !********************************************************************** subroutine InitAll - ! System intialization + ! System initialization call InitMPI() call GaussQuad_Init call ReadConfig('Input/tube.in') diff --git a/examples/case/Makefile b/examples/case/Makefile index d9307f7..ac1458d 100755 --- a/examples/case/Makefile +++ b/examples/case/Makefile @@ -1,22 +1 @@ -include ../../Makefile.in - -LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(LAPACK_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC) -EXECUTABLES = tube initcond - -all : $(EXECUTABLES) - -lib : $(wildcard $(WORK_DIR)/common/*.F90) - make -C $(WORK_DIR)/common - -$(EXECUTABLES) : % : %.o $(wildcard $(COMMON_DIR)/*.F90) - make lib - $(FC) $(LDFLAGS) -o $@ $< $(LIB) - -clean : - -rm -f $(EXECUTABLES) *.o *.mod *.a core - -# Dependency -.depend : $(wildcard *.F90) - $(MAKEDEPEND_BIN) $^ > .depend - -include .depend +include ../Makefile.mkl diff --git a/examples/case/initcond.F90 b/examples/case/initcond.F90 index 245bd48..82ead4c 100644 --- a/examples/case/initcond.F90 +++ b/examples/case/initcond.F90 @@ -105,7 +105,7 @@ program InitCond Nt0 = 0; time = 0. vBkg(1:2) = 0.; vBkg(3) = 8. - ! Write intial conditions + ! Write initial conditions if (nrbc > 0) then write (fn, FMT=fn_FMT) 'D/', 'x', 0, '.dat' call WriteManyRBCs(fn, nrbc, rbcs) diff --git a/examples/case/run.sh b/examples/case/run.sh new file mode 100644 index 0000000..6194ad3 --- /dev/null +++ b/examples/case/run.sh @@ -0,0 +1,34 @@ +#!/bin/bash + +#SBATCH --account=gts-sbryngelson3 +#SBATCH -N4 --ntasks-per-node=24 +#SBATCH --mem-per-cpu=2G +#SBATCH -t1:00:00 +#SBATCH -q embers +#SBATCH --mail-user=smanasreh6@gatech.edu +#SBATCH --mail-type=BEGIN,END,FAIL +#SBATCH -o "./run_logs/comeon.log" + +cd $SLURM_SUBMIT_DIR + +cd D +rm -rf *x* +rm -rf r* +rm -rf w* +cd ../ + +ml gcc mvapich2 mkl netcdf-c netcdf-cxx netcdf-fortran fftw + +cd ../../common +make clean +make .depend +make + +cd ../examples/case +make clean +make .depend +make +srun -n 1 ./initcond +srun ./tube + +# sbatch run.sh \ No newline at end of file diff --git a/examples/case/tube.F90 b/examples/case/tube.F90 index c68db12..73f0961 100755 --- a/examples/case/tube.F90 +++ b/examples/case/tube.F90 @@ -12,12 +12,13 @@ program cells_in_tube use ModBasicMath use ModQuadRule +#include "petsc/finclude/petsc.h" + use petsc + implicit none integer :: cutoff character(CHRLEN) :: fn -#include "../../petsc_include.h" - call InitAll call TimeInt_Euler @@ -31,7 +32,7 @@ program cells_in_tube !********************************************************************** subroutine InitAll - ! System intialization + ! System initialization call InitMPI() call GaussQuad_Init call ReadConfig('Input/tube.in') diff --git a/examples/case_diff_celltypes/Makefile b/examples/case_diff_celltypes/Makefile index d9307f7..c68a502 100755 --- a/examples/case_diff_celltypes/Makefile +++ b/examples/case_diff_celltypes/Makefile @@ -1,22 +1 @@ -include ../../Makefile.in - -LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(LAPACK_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC) -EXECUTABLES = tube initcond - -all : $(EXECUTABLES) - -lib : $(wildcard $(WORK_DIR)/common/*.F90) - make -C $(WORK_DIR)/common - -$(EXECUTABLES) : % : %.o $(wildcard $(COMMON_DIR)/*.F90) - make lib - $(FC) $(LDFLAGS) -o $@ $< $(LIB) - -clean : - -rm -f $(EXECUTABLES) *.o *.mod *.a core - -# Dependency -.depend : $(wildcard *.F90) - $(MAKEDEPEND_BIN) $^ > .depend - -include .depend +include ../Makefile.mkl \ No newline at end of file diff --git a/examples/case_diff_celltypes/initcond.F90 b/examples/case_diff_celltypes/initcond.F90 index df13da6..a52a866 100644 --- a/examples/case_diff_celltypes/initcond.F90 +++ b/examples/case_diff_celltypes/initcond.F90 @@ -116,7 +116,7 @@ program InitCond Nt0 = 0; time = 0. vBkg(1:2) = 0.; vBkg(3) = 8. - ! Write intial conditions + ! Write initial conditions if (nrbc > 0) then write (fn, FMT=fn_FMT) 'D/', 'x', 0, '.dat' call WriteManyRBCs(fn, nrbc, rbcs) diff --git a/examples/case_diff_celltypes/tube.F90 b/examples/case_diff_celltypes/tube.F90 index 61160da..3f492ca 100755 --- a/examples/case_diff_celltypes/tube.F90 +++ b/examples/case_diff_celltypes/tube.F90 @@ -12,12 +12,13 @@ program cells_in_tube use ModBasicMath use ModQuadRule +#include "petsc/finclude/petsc.h" + use petsc + implicit none integer :: cutoff character(CHRLEN) :: fn -#include "../../petsc_include.h" - call InitAll call TimeInt_Euler @@ -31,7 +32,7 @@ program cells_in_tube !********************************************************************** subroutine InitAll - ! System intialization + ! System initialization call InitMPI() call GaussQuad_Init call ReadConfig('Input/tube.in') diff --git a/examples/field_visual/README.md b/examples/field_visual/README.md index dd26c44..b490093 100644 --- a/examples/field_visual/README.md +++ b/examples/field_visual/README.md @@ -16,7 +16,7 @@ Then, modify the `./Input/tube.in` configuration file to match the restart file' ### Output File -After runnning the simulation, you will receive an output file `./D/field.csv`. This file contains the following comma-separated columns: +After running the simulation, you will receive an output file `./D/field.csv`. This file contains the following comma-separated columns: X | Y | Z | Vx | Vy | Vz diff --git a/examples/field_visual/field.F90 b/examples/field_visual/field.F90 index f050df6..99ecaeb 100644 --- a/examples/field_visual/field.F90 +++ b/examples/field_visual/field.F90 @@ -173,7 +173,7 @@ END SUBROUTINE leukVel !********************************************************************** subroutine InitAll - ! System intialization + ! System initialization call InitMPI() call GaussQuad_Init call ReadConfig('Input/tube.in') diff --git a/examples/field_visual/traction.F90 b/examples/field_visual/traction.F90 index 020ddf5..09deb8d 100644 --- a/examples/field_visual/traction.F90 +++ b/examples/field_visual/traction.F90 @@ -54,7 +54,7 @@ end subroutine GetTraction !********************************************************************** subroutine InitAll - ! System intialization + ! System initialization call InitMPI() call GaussQuad_Init call ReadConfig('Input/tube.in') diff --git a/examples/randomized_case/Makefile b/examples/randomized_case/Makefile index 149406a..c68a502 100755 --- a/examples/randomized_case/Makefile +++ b/examples/randomized_case/Makefile @@ -1,22 +1 @@ -include ../../Makefile.in - -LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(LAPACK_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC) -EXECUTABLES = initcond tube - -all : $(EXECUTABLES) - -lib : $(wildcard $(WORK_DIR)/common/*.F90) - make -C $(WORK_DIR)/common - -$(EXECUTABLES) : % : %.o $(wildcard $(COMMON_DIR)/*.F90) - make lib - $(FC) $(LDFLAGS) -o $@ $< $(LIB) - -clean : - -rm -f $(EXECUTABLES) *.o *.mod *.a core - -# Dependency -.depend : $(wildcard *.F90) - $(MAKEDEPEND_BIN) $^ > .depend - -include .depend +include ../Makefile.mkl \ No newline at end of file diff --git a/examples/randomized_case/initcond.F90 b/examples/randomized_case/initcond.F90 index 7d53b2b..a8c3f07 100644 --- a/examples/randomized_case/initcond.F90 +++ b/examples/randomized_case/initcond.F90 @@ -244,7 +244,7 @@ subroutine rotate_cell(cell) zv(2) = v2*2*sqrt(1 - vsq) zv(3) = 1 - (2*vsq) - !generate rotation matrix + !generate rotation matrix R rotmat = RotateMatrix(zv) !(/ 0., 1., 0./)) !rotate cell with rotmat @@ -280,7 +280,7 @@ logical function check_wall_collision(cell) sp = walls(i2)%x(j2, :) !just do a distance check - if (VecNorm(sp - cp) .le. threshold) then + if (NORM2(sp - cp) .le. threshold) then check_wall_collision = .true. return end if @@ -307,7 +307,7 @@ logical function check_cell_collision(cell1, cell2) check_cell_collision = .false. !first check the centers; if distance > 4 then we are sure that they don't collide - if (VecNorm(cell1%xc - cell2%xc) .ge. 4) return + if (NORM2(cell1%xc - cell2%xc) .ge. 4) return !each point in cell1 do i = 1, cell1%nlat diff --git a/examples/randomized_case/tube.F90 b/examples/randomized_case/tube.F90 index 9e2f7c3..ccbc49c 100755 --- a/examples/randomized_case/tube.F90 +++ b/examples/randomized_case/tube.F90 @@ -12,12 +12,13 @@ program cells_in_tube use ModBasicMath use ModQuadRule +#include "petsc/finclude/petsc.h" + use petsc + implicit none integer :: cutoff character(CHRLEN) :: fn -#include "../../petsc_include.h" - call InitAll call TimeInt_Euler @@ -31,7 +32,7 @@ program cells_in_tube !********************************************************************** subroutine InitAll - ! System intialization + ! System initialization call InitMPI() call GaussQuad_Init call ReadConfig('Input/tube.in') diff --git a/install/format.sh b/install/format.sh new file mode 100644 index 0000000..8c4fa9a --- /dev/null +++ b/install/format.sh @@ -0,0 +1,10 @@ +#!/bin/bash + +echo "Formatting example cases and common directory in $PWD" + +pip3 install --upgrade fprettify + +fprettify ./examples -r --indent 2 +fprettify ./common -r --indent 2 + +echo "Done formatting!" \ No newline at end of file diff --git a/install/install-with-lapack.sh b/install/install-with-lapack.sh new file mode 100644 index 0000000..b781408 --- /dev/null +++ b/install/install-with-lapack.sh @@ -0,0 +1,48 @@ +#!/bin/bash + +# salloc a node before you run this because petsc configure uses srun + +# replace with equivalent modules on your cluster +# if module is not available, follow readme.md to manually install +ml gcc mvapich2 python/3.9.12-rkxvr6 netcdf-c netcdf-cxx netcdf-fortran fftw + +# create packages directory +mkdir packages +cd packages + +# build and install lapack and blas +wget https://github.com/Reference-LAPACK/lapack/archive/refs/tags/v3.11.tar.gz +tar -xvf v3.11.tar.gz +cd lapack-3.11 +cp ../../install/scripts/make.inc ./ +make + +# build and install petsc 3.19.6 in packages directory +wget https://ftp.mcs.anl.gov/pub/petsc/petsc-3.19.tar.gz +tar -xvf petsc-3.19.tar.gz + +pip3 install --user configure + +cp ../install/scripts/petsc_configure.py ./petsc-3.19.6 +cd petsc-3.19.6 +python3 petsc_configure.py --blas-lapack --dryrun +python3 petsc_configure.py --blas-lapack + +make PETSC_DIR=`pwd` PETSC_ARCH=petsc_configure all +make PETSC_DIR=`pwd` PETSC_ARCH=petsc_configure check + +# build and install spherepack +cd .. +git clone https://github.com/comp-physics/spherepack3.2.git +cd spherepack3.2 +make + +# build and install makedepf90 +cd .. +git clone https://github.com/comp-physics/makedepf90.git +cd makedepf90 +# it works without setting prefix +make +make install + +echo "Done installing!" \ No newline at end of file diff --git a/install/install-with-mkl.sh b/install/install-with-mkl.sh new file mode 100644 index 0000000..ef504a2 --- /dev/null +++ b/install/install-with-mkl.sh @@ -0,0 +1,40 @@ +#!/bin/bash + +# salloc a node before you run this because petsc configure uses srun + +# replace with equivalent modules on your cluster +# if module is not available, follow readme.md to manually install +ml gcc mvapich2 mkl python/3.9.12-rkxvr6 netcdf-c netcdf-cxx netcdf-fortran fftw + +# building and installing petsc 3.19.6 in packages directory +mkdir packages +cd packages + +wget https://ftp.mcs.anl.gov/pub/petsc/petsc-3.19.tar.gz +tar -xvf petsc-3.19.tar.gz + +pip3 install --user configure + +cp ../install/scripts/petsc_configure.py ./petsc-3.19.6 +cd petsc-3.19.6 +python3 petsc_configure.py --mkl-only --dryrun +python3 petsc_configure.py --mkl-only + +make PETSC_DIR=`pwd` PETSC_ARCH=petsc_configure all +make PETSC_DIR=`pwd` PETSC_ARCH=petsc_configure check + +# build and install spherepack +cd .. +git clone https://github.com/comp-physics/spherepack3.2.git +cd spherepack3.2 +make + +# build and install makedepf90 +cd .. +git clone https://github.com/comp-physics/makedepf90.git +cd makedepf90 +# it works without setting prefix +make +make install + +echo "Done installing!" \ No newline at end of file diff --git a/install/readme.md b/install/readme.md index f20a2fd..41c330d 100644 --- a/install/readme.md +++ b/install/readme.md @@ -21,6 +21,12 @@ You will need `gcc`, `gfortran`, and a suitable MPI wrapper like `mvapich` (or t * To check for gfortran and gcc, see if `which gfortran` and `which gcc` return a path * Similarly, to see if you can run MPI commands for later, see if `which mpicc` or `which mpif90` return a path +You will also need python3 for the PETSc install and pip modules. On Phoenix, you can module load it via `ml python/3.9.12-rkxvr6`. On Phoenix, you may also need to add the `~/.local/bin` directory to your PATH by adding this line to your `~/.bashrc`: + +```shell +export PATH="$PATH:$HOME/.local/bin" +``` + ## Build libraries ### MKL @@ -29,22 +35,10 @@ You will need `gcc`, `gfortran`, and a suitable MPI wrapper like `mvapich` (or t * This will automatically set the `MKL_ROOT` environment variable necessary for `Makefile.in` * You can check environment variables via `module show mkl` * If your cluster does not have mkl, it's available for download and install [here](https://www.intel.com/content/www/us/en/developer/tools/oneapi/base-toolkit-download.html?operatingsystem=linux&distributions=offline). -* For a manual mkl install, we recommend installing inside a packages directory. Creating an `MKL_ROOT` variable in `Makefile.in` may be necessary too. It should be set to the path of the `/mkl` directory. +* If MKL is not available on your system, follow the LAPACK and BLAS step, and skip this step. * Note that `MKL_LIB` options in Makefile.in may need to be changed depending on the version of mkl, but the [mkl link line advisor](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-link-line-advisor.html#gs.9hbhed) should provide the correct link options. -### BLAS - -* Make a packages directory inside the RBC3D directory: `mkdir packages` -* Move back into the packages: `cd packages` -* Download the latest BLAS (at time of writing `3.11.0`): `wget http://www.netlib.org/blas/blas-3.11.0.tgz` -* Unpack it: `tar -xvf blas-3.11.0.tgz` -* `cd BLAS-3.11.0` -* Modify `make.inc` line 18 as `FC = mpif90` -* Execute `make`, which will create the library file `blas_LINUX.a` -* Later, You will need the absolute path of `blas_LINUX.a` to configure `petsc-lite` - * In my case this is `/storage/coda1/p-sbryngelson3/0/sbryngelson3/RBC3D/packages/BLAS-3.11.0/blas_LINUX.a` - -### LAPACK +### LAPACK and BLAS * Move back into packages: `cd RBC3D/packages` * Download the latest lapack (at time of writing `3.11`): `wget https://github.com/Reference-LAPACK/lapack/archive/refs/tags/v3.11.tar.gz` @@ -55,37 +49,31 @@ You will need `gcc`, `gfortran`, and a suitable MPI wrapper like `mvapich` (or t * `FC = mpif90` (line 20) * `mv make.inc.example make.inc` * Build (this takes a few minutes): `make` -* Later, You will need the absolute path of `liblapack.a` to configure `petsc-lite` +* Later, You will need the absolute path of `liblapack.a` and `librefblas.a` to configure `petsc-lite` * In my case this is `/storage/coda1/p-sbryngelson3/0/sbryngelson3/RBC3D/packages/lapack-3.11/liblapack.a` + * and `/storage/coda1/p-sbryngelson3/0/sbryngelson3/RBC3D/packages/lapack-3.11/librefblas.a` +* Note that if you're choosing to install LAPACK/BLAS instead of MKL, you'll need to include `Makefile.lapack` in `examples/case/Makefile` instead of `Makefile.mkl` when you run the example case later. -### Valgrind -* Easiest if this is already available or can be loaded. -* PACE Phoenix has this available as a module: `module load valgrind` -* `module show valgrind` or `which valgrind` can tell you where the library is. -* At time of writing, it is here: `/usr/local/pace-apps/manual/packages/valgrind/3.19.0/gcc-4.8.5` - * You will need this path to build PETSc-lite - -### PETSc-lite +### PETSc -* This depends on Valgrind, LAPACK, BLAS from above, don't attempt until those steps are finished +* This depends on MKL or LAPACK/BLAS from above, don't attempt until one of those steps is finished * Move back to `RBC3D/packages` -* Copy `petsc-lite-3.0.0-p3.tar.gz` from the libs directory into the packages directory: `cp ../libs/petsc-lite-3.0.0-p3.tar.gz ./` -* Unpack `petsc-lite`: `tar -xvf petsc-lite-3.0.0-p3.tar.gz` -* Set up your environment via the environment variables - * Get the absolute path of the unpacked petsc via - * `cd petsc-3.0.0-p3` - * `pwd`: In my case: `/storage/home/hcoda1/6/sbryngelson3/p-sbryngelson3-0/RBC3D/packages/petsc-3.0.0-p3` - * If you are using something other than bash, look up how to set environment variables for it, otherwise: - * Execute (notice this is the path from above): `export PETSC_DIR=/storage/home/hcoda1/6/sbryngelson3/p-sbryngelson3-0/RBC3D/packages/petsc-3.0.0-p3` - * Execute: `export PETSC_ARCH=linux-gnu-c-opt` -* Ascend up a directory and create a new build directory like `mkdir RBC3D/packages/mypetsc` then cd back into `petsc-3.0.0-p3` -* Configure via something like this, using your own absolute paths (for blas, lapack, valgrind, and mypetsc), and notice the `--with-mpiexec=srun` line where you should replace `srun` with what is relevant for your system (`srun` if available, `mpirun` or `mpiexec` are two other options): -``` -./configure --with-cc=mpicc --with-fc=mpif90 --with-debugging=0 COPTFLAGS='-O3 -march=native -mtune=native' CXXOPTFLAGS='-O3 -march=native -mtune=native' FOPTFLAGS='-O3 -march=native -mtune=native' --with-blas-lib=/storage/coda1/p-sbryngelson3/0/sbryngelson3/RBC3D/packages/BLAS-3.11.0/blas_LINUX.a --with-lapack-lib=/storage/coda1/p-sbryngelson3/0/sbryngelson3/RBC3D/packages/lapack-3.11/liblapack.a --with-valgrind-dir=/usr/local/pace-apps/manual/packages/valgrind/3.19.0/gcc-4.8.5 --prefix=/storage/coda1/p-sbryngelson3/0/sbryngelson3/RBC3D/packages/mypetsc --with-shared=0 --with-mpiexec=srun --with-x11=0 --with-x=0 --with-windows-graphics=0 -``` -* Build: `make` -* Install: `make install` +* Download PETSc: wget https://ftp.mcs.anl.gov/pub/petsc/petsc-3.19.tar.gz` +* Unpack `petsc-3.19`: `tar -xvf petsc-3.19.tar.gz` +* Copy configure script into petsc directory: `cp ../install/py_scripts/petsc_configure.py ./petsc-3.19.6` +* Descend into the directory: `cd petsc-3.19` +* Install pip configure: `pip3 install --user configure` +* If you installed MKL, run: `python3 petsc_configure.py --mkl-only` +* If you installed LAPACK/BLAS instead, run: `python3 petsc_configure.py --blas-lapack` +* Note the `--dryrun` option will show you the PETSc configure options. + +* Build and Test: +```shell +make PETSC_DIR=`pwd` PETSC_ARCH=petsc_configure all +make PETSC_DIR=`pwd` PETSC_ARCH=petsc_configure check +``` +* `make check` may return errors, but you can ignore these if the tests completed. ### FFTW diff --git a/install/scripts/make.inc b/install/scripts/make.inc new file mode 100644 index 0000000..6db4b49 --- /dev/null +++ b/install/scripts/make.inc @@ -0,0 +1,86 @@ +#################################################################### +# LAPACK make include file. # +#################################################################### + +SHELL = /bin/sh + +# CC is the C compiler, normally invoked with options CFLAGS. +# +CC = mpicc +CFLAGS = -O3 + +# Modify the FC and FFLAGS definitions to the desired compiler +# and desired compiler options for your machine. NOOPT refers to +# the compiler options desired when NO OPTIMIZATION is selected. +# +# Note: During a regular execution, LAPACK might create NaN and Inf +# and handle these quantities appropriately. As a consequence, one +# should not compile LAPACK with flags such as -ffpe-trap=overflow. +# +FC = mpif90 +FFLAGS = -O2 -frecursive +FFLAGS_DRV = $(FFLAGS) +FFLAGS_NOOPT = -O0 -frecursive + +# Define LDFLAGS to the desired linker options for your machine. +# +LDFLAGS = + +# The archiver and the flag(s) to use when building an archive +# (library). If your system has no ranlib, set RANLIB = echo. +# +AR = ar +ARFLAGS = cr +RANLIB = ranlib + +# Timer for the SECOND and DSECND routines +# +# Default: SECOND and DSECND will use a call to the +# EXTERNAL FUNCTION ETIME +#TIMER = EXT_ETIME +# For RS6K: SECOND and DSECND will use a call to the +# EXTERNAL FUNCTION ETIME_ +#TIMER = EXT_ETIME_ +# For gfortran compiler: SECOND and DSECND will use a call to the +# INTERNAL FUNCTION ETIME +TIMER = INT_ETIME +# If your Fortran compiler does not provide etime (like Nag Fortran +# Compiler, etc...) SECOND and DSECND will use a call to the +# INTERNAL FUNCTION CPU_TIME +#TIMER = INT_CPU_TIME +# If none of these work, you can use the NONE value. +# In that case, SECOND and DSECND will always return 0. +#TIMER = NONE + +# Uncomment the following line to include deprecated routines in +# the LAPACK library. +# +#BUILD_DEPRECATED = Yes + +# LAPACKE has the interface to some routines from tmglib. +# If LAPACKE_WITH_TMG is defined, add those routines to LAPACKE. +# +#LAPACKE_WITH_TMG = Yes + +# Location of the extended-precision BLAS (XBLAS) Fortran library +# used for building and testing extended-precision routines. The +# relevant routines will be compiled and XBLAS will be linked only +# if USEXBLAS is defined. +# +#USEXBLAS = Yes +#XBLASLIB = -lxblas + +# The location of the libraries to which you will link. (The +# machine-specific, optimized BLAS library should be used whenever +# possible.) +# +BLASLIB = $(TOPSRCDIR)/librefblas.a +CBLASLIB = $(TOPSRCDIR)/libcblas.a +LAPACKLIB = $(TOPSRCDIR)/liblapack.a +TMGLIB = $(TOPSRCDIR)/libtmglib.a +LAPACKELIB = $(TOPSRCDIR)/liblapacke.a + +# DOCUMENTATION DIRECTORY +# If you generate html pages (make html), documentation will be placed in $(DOCSDIR)/explore-html +# If you generate man pages (make man), documentation will be placed in $(DOCSDIR)/man +DOCSDIR = $(TOPSRCDIR)/DOCS diff --git a/install/scripts/petsc_configure.py b/install/scripts/petsc_configure.py new file mode 100644 index 0000000..c376ffd --- /dev/null +++ b/install/scripts/petsc_configure.py @@ -0,0 +1,144 @@ +#!/usr/bin/env python +""" Configuration helper for PETSc + +Execute from PETSC_DIR + +Note that older versions of PETSc require Python 2 + +Run with -h to see arguments + +Adds a thin extra layer around PETSc's configuration script, +to help combine commonly-used combinations of configuration options and +construct meaningful PETSC_ARCH values. + +Proceeds by collecting the set of passed arguments and processing them +before sending them to PETSc's configure script. +This logic will likely be somewhat brittle. Always do a sanity check +and look at the options that are actually being sent. This script should +be simple enough to figure out what's going on. +""" + +from __future__ import print_function +import sys +import os +import argparse +import re + + +def main(): + """ Main script logic """ + args, options_in = get_args() + options = process_args(options_in, args) + petsc_configure(options, args) + + +def get_args(): + """ Retrieve custom arguments and remaining arguments""" + parser = argparse.ArgumentParser( + description='Compute arguments to pass to PETSc\'s configure script') + parser.add_argument( + '--dryrun', + action="store_true", + help="don't actually configure") + parser.add_argument( + '--mkl-only', + action="store_true", + help="Configure with MKL") + parser.add_argument( + '--blas-lapack', + action="store_true", + help="Configure with blas/lapack") + args, unknown = parser.parse_known_args() + return args, unknown + +def process_args(options_in, args): + """ Main logic to create a set of options for PETSc's configure script, + along with a corresponding PETSC_ARCH string, if required """ + + if args.mkl_only: + return options_for_mkl() + + if args.blas_lapack: + return options_for_blaslapack() + + +def options_for_mkl(): + """ Return a custom set of arguments to simply download and build MPICH """ + + MKL_DIR = os.getenv('MKLROOT') + + options = [] + options.append('--with-cc=mpicc') + options.append('--with-cxx=mpicxx') + options.append('--with-fc=mpif90') + options.append('--with-fortran-datatypes') + options.append('--with-debugging=0') + options.append("--COPTFLAGS=-g -O3 -march=native") + options.append("--CXXOPTFLAGS=-g -O3 -march=native") + options.append("--FOPTFLAGS=-g -O3 -march=native") + options.append("--with-blaslapack-dir=" + MKL_DIR) + options.append("--with-mpiexec=srun") + # this is an option from the original configure but i think this is unnecessary + # options.append("--with-shared-libraries=0") + options.append("--with-x11=0") + options.append("--with-x=0") + options.append("--with-windows-graphics=0") + return options + +def options_for_blaslapack(): + """ Return a custom set of arguments to simply download and build MPICH """ + + PACKAGES_DIR = os.path.dirname(os.getcwd()) + LAPACK_LIB = PACKAGES_DIR + "/lapack-3.11/liblapack.a" + BLAS_LIB = PACKAGES_DIR + "/lapack-3.11/librefblas.a" + + print("PACKAGES_DIR: " + PACKAGES_DIR) + print("BLAS_LIB: " + BLAS_LIB) + print("LAPACK_LIB: " + LAPACK_LIB) + + options = [] + options.append('--with-cc=mpicc') + options.append('--with-cxx=mpicxx') + options.append('--with-fc=mpif90') + options.append('--with-fortran-datatypes') + options.append('--with-debugging=0') + options.append("--COPTFLAGS=-g -O3 -march=native") + options.append("--CXXOPTFLAGS=-g -O3 -march=native") + options.append("--FOPTFLAGS=-g -O3 -march=native") + options.append("--with-blas-lib=" + BLAS_LIB) + options.append("--with-lapack-lib=" + LAPACK_LIB) + options.append("--with-mpiexec=srun") + # configure throws an error if the option below isn't on? + options.append("--with-shared-libraries=0") + options.append("--with-x11=0") + options.append("--with-x=0") + options.append("--with-windows-graphics=0") + return options + +def petsc_configure(options, args): + """ Standard PETSc configuration script logic (from config/examples) """ + if args.dryrun: + print("Dry Run. Would configure with these options:") + print("\n".join(options)) + else: + sys.path.insert(0, os.path.abspath('config')) + # print(sys.version) + try: + import configure + except ImportError: + print( + 'PETSc configure module not found. Make sure you are executing from PETSC_DIR' + ) + sys.exit(1) + print('Configuring with these options:') + print("\n".join(options)) + + # Since petsc_configure looks directly at sys.argv, remove and replace arguments + argv_temp = sys.argv + sys.argv = sys.argv[:1] + configure.petsc_configure(options) + sys.argv = argv_temp + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/libs/petsc-lite-3.0.0-p3.tar.gz b/libs/petsc-lite-3.0.0-p3.tar.gz deleted file mode 100644 index fadce42..0000000 Binary files a/libs/petsc-lite-3.0.0-p3.tar.gz and /dev/null differ diff --git a/petsc_include.h b/petsc_include.h deleted file mode 100755 index 601698b..0000000 --- a/petsc_include.h +++ /dev/null @@ -1,6 +0,0 @@ -#define PETSC_AVOID_MPIF_H -#include "finclude/petsc.h" -#include "finclude/petscvec.h" -#include "finclude/petscmat.h" -#include "finclude/petscksp.h" -#include "finclude/petscpc.h" diff --git a/rbc.sh b/rbc.sh new file mode 100644 index 0000000..661d229 --- /dev/null +++ b/rbc.sh @@ -0,0 +1,33 @@ +#!/bin/bash + +# change this to the correct module load for python3 on your cluster +ml python/3.9.12-rkxvr6 + +if [[ ":$PATH:" != *":$HOME/.local/bin:"* ]]; then + export PATH="$PATH:$HOME/.local/bin" +fi + +# check python3 works +command -v python3 > /dev/null 2>&1 +if (($?)); then + echo "[rbc.sh] Error: Couldn't find Python3. Please ensure it is discoverable." + exit 1 +fi + +python3 -c 'print("")' > /dev/null +if (($?)); then + echo "[rbc.sh] Error: Python3 is present but can't execute a simple program. Please ensure that python3 is working." + exit 1 +fi + +if [ "$1" == 'format' ]; then + . "$(pwd)/install/format.sh" $@; exit +fi + +if [ "$1" == 'install-with-mkl' ]; then + . "$(pwd)/install/install-with-mkl.sh" $@; exit +fi + +if [ "$1" == 'install-with-lapack' ]; then + . "$(pwd)/install/install-with-lapack.sh" $@; exit +fi \ No newline at end of file