Merge pull request #73 from Treecodes/develop

Merge in RBS kernels

examples/random_cube_reproducible.c

@@ -171,16 +171,16 @@ int main(int argc, char **argv)
                 double xx, yy, zz;
                 Point_Plummer_Octant(plummer_R , &xx, &yy, &zz);
 
-		for (int ii = 0; ii < 2; ++ii) {
-		    for (int jj = 0; jj < 2; ++jj) {
-			for (int kk = 0; kk < 2; ++kk) {
-		    	    int index = (N/8) * (ii*4 + jj*2 + kk) + i;
-		    	    mySources.x[index] = xx * pow(-1, ii);
-		            mySources.y[index] = yy * pow(-1, jj);
-		            mySources.z[index] = zz * pow(-1, kk);
-                	}
-		    }
-		}
+		        for (int ii = 0; ii < 2; ++ii) {
+		            for (int jj = 0; jj < 2; ++jj) {
+			            for (int kk = 0; kk < 2; ++kk) {
+		    	            int index = (N/8) * (ii*4 + jj*2 + kk) + i;
+		    	            mySources.x[index] = xx * pow(-1, ii);
+		                    mySources.y[index] = yy * pow(-1, jj);
+		                    mySources.z[index] = zz * pow(-1, kk);
+              	        }
+		            }
+		        }
             }
         }
 
@@ -208,6 +208,46 @@ int main(int argc, char **argv)
             }
         }
 
+    } else if (distribution == SLAB_1) {
+
+        for (int j = 0; j < rank+1; ++j) { // Cycle to generate same particle no matter num ranks
+            for (int i = 0; i < N; ++i) {
+                mySources.x[i] = ((double)random()/(double)(RAND_MAX)) * 10. - 5;
+                mySources.y[i] = ((double)random()/(double)(RAND_MAX)) * 10. - 5.;
+                mySources.z[i] = ((double)random()/(double)(RAND_MAX)) * 1.;
+                mySources.q[i] = ((double)random()/(double)(RAND_MAX)) * 2. - 1.;
+                mySources.w[i] = 1.0;
+                mySources.myGlobalIDs[i] = (ZOLTAN_ID_TYPE)(rank*N + i);
+                mySources.b[i] = 1.0; // dummy weighting scheme
+            }
+        }
+
+    } else if (distribution == SLAB_2) {
+
+        for (int j = 0; j < rank+1; ++j) { // Cycle to generate same particle no matter num ranks
+            for (int i = 0; i < N; ++i) {
+                mySources.x[i] = ((double)random()/(double)(RAND_MAX)) * 10 - 5.;
+                mySources.y[i] = ((double)random()/(double)(RAND_MAX)) * 1.;
+                mySources.z[i] = ((double)random()/(double)(RAND_MAX)) * 1.;
+                mySources.q[i] = ((double)random()/(double)(RAND_MAX)) * 2. - 1.;
+                mySources.w[i] = 1.0;
+                mySources.myGlobalIDs[i] = (ZOLTAN_ID_TYPE)(rank*N + i);
+                mySources.b[i] = 1.0; // dummy weighting scheme
+            }
+        }
+
+    } else if (distribution == SPHERICAL_SHELL) {
+
+        for (int j = 0; j < rank+1; ++j) { //Cycle to generate same particle no matter num ranks
+            for (int i = 0; i < N; ++i) {
+                Point_Spherical_Shell(1., &mySources.x[i], &mySources.y[i], &mySources.z[i]);
+                mySources.q[i] = ((double)random()/(double)(RAND_MAX)) * 2. - 1.;
+                mySources.w[i] = 1.0;
+                mySources.myGlobalIDs[i] = (ZOLTAN_ID_TYPE)(rank*N + i);
+                mySources.b[i] = 1.0;
+            }
+        }
+
     } else {
         printf("[random cube example] ERROR! Distribution %d undefined in this "
                 "context.  Exiting.\n", distribution);

examples/support_fns.c

@@ -265,6 +265,15 @@ void Params_Parse(FILE *fp, struct RunParams **run_params, int *N, int *M, int *
     } else if (strcasecmp(distribution_string, "PLUMMER_SYMMETRIC") == 0) {
         *distribution = PLUMMER_SYMMETRIC;
 
+    } else if (strcasecmp(distribution_string, "SLAB_1") == 0) {
+        *distribution = SLAB_1;
+
+    } else if (strcasecmp(distribution_string, "SLAB_2") == 0) {
+        *distribution = SLAB_2;
+
+    } else if (strcasecmp(distribution_string, "SPHERICAL_SHELL") == 0) {
+        *distribution = SPHERICAL_SHELL;
+
     } else {
         if (rank == 0) {
             printf("[random cube example] ERROR! Undefined distribution token \"%s\". Exiting.\n",
@@ -443,6 +452,24 @@ void Point_Exponential(double *x, double *y, double *z)
 }
 
 
+/*----------------------------------------------------------------------------*/
+void Point_Spherical_Shell(double R, double *x, double *y, double *z)
+{
+    double u = 2. * M_PI * (double)random()/(1.+ (double)(RAND_MAX));
+    double v = 2. * (double)random()/(1.+ (double)(RAND_MAX)) - 1.;
+
+    *x = sqrt(1. - v * v) * cos(u);
+    *y = sqrt(1. - v * v) * sin(u);
+    *z = v;
+
+    *x *= R;
+    *y *= R;
+    *z *= R;
+
+    return;
+}
+
+
 
 /*----------------------------------------------------------------------------*/
 void Timing_Calculate(double time_run_glob[3][4], double time_tree_glob[3][13], double time_direct_glob[3][4],

examples/support_fns.h

@@ -13,7 +13,10 @@ typedef enum DISTRIBUTION
     GAUSSIAN,
     EXPONENTIAL,
     PLUMMER,
-    PLUMMER_SYMMETRIC
+    PLUMMER_SYMMETRIC,
+    SLAB_1,
+    SLAB_2,
+    SPHERICAL_SHELL
 } DISTRIBUTION;
 
 typedef enum PARTITION
@@ -40,6 +43,8 @@ void Point_Gaussian(double *x, double *y, double *z);
 
 void Point_Exponential(double *x, double *y, double *z);
 
+void Point_Spherical_Shell(double R, double *x, double *y, double *z);
+
 
 void Timing_Calculate(double time_run_glob[3][4], double time_tree_glob[3][13], double time_direct_glob[3][4],
                 double time_run[4], double time_tree[13], double time_direct[4]);

interfaces/fortran/BaryTreeInterface.fh

@@ -0,0 +1,39 @@
+ENUM, BIND(C)
+    ENUMERATOR :: NO_KERNEL = 0, COULOMB, YUKAWA, &
+                  REGULARIZED_COULOMB, REGULARIZED_YUKAWA, &
+                  ATANF, TCF, DCF, SIN_OVER_R, MQ, RBS_U, RBS_V, USER
+
+    ENUMERATOR :: NO_SINGULARITY = 0, SKIPPING, SUBTRACTION
+
+    ENUMERATOR :: NO_APPROX = 0, LAGRANGE, HERMITE
+
+    ENUMERATOR :: NO_COMPUTE_TYPE = 0, PARTICLE_CLUSTER, &
+                  CLUSTER_PARTICLE, CLUSTER_CLUSTER
+END ENUM
+
+INTERFACE
+    SUBROUTINE BaryTreeInterface(num_targets, num_sources, &
+                    target_x, target_y, target_z, target_q, &
+                    source_x, source_y, source_z, source_q, source_w, &
+                    potential, kernel, num_kernel_params, kernel_params, &
+                    singularity, approximation, compute_type, &
+                    theta, degree, max_source_leaf, max_target_leaf, &
+                    size_check, beta, verbosity) &
+               BIND(C, NAME='BaryTreeInterface')
+
+        USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT, C_DOUBLE, C_PTR
+        IMPLICIT NONE
+
+        INTEGER(KIND=C_INT), VALUE, INTENT(IN) :: num_targets, num_sources
+        TYPE(C_PTR), VALUE, INTENT(IN) :: target_x, target_y, target_z, target_q
+        TYPE(C_PTR), VALUE, INTENT(IN) :: source_x, source_y, source_z, source_q, source_w
+        TYPE(C_PTR), VALUE, INTENT(IN) :: potential
+
+        INTEGER(KIND=C_INT), VALUE, INTENT(IN) :: kernel, num_kernel_params, &
+                            singularity, approximation, compute_type, degree, &
+                            max_source_leaf, max_target_leaf, verbosity
+        TYPE(C_PTR), VALUE, INTENT(IN) :: kernel_params
+        REAL(KIND=C_DOUBLE), VALUE, INTENT(IN) :: theta, size_check, beta
+
+    END SUBROUTINE BaryTreeInterface
+END INTERFACE

interfaces/fortran/example_interface.f03

@@ -0,0 +1,91 @@
+PROGRAM BaryTree_Fortran_Example
+
+    USE MPI
+    USE ISO_C_BINDING, ONLY: C_LOC
+    IMPLICIT NONE
+
+    INCLUDE "BaryTreeInterface.fh"
+
+    INTEGER :: rank, num_proc, ierr
+
+    INTEGER :: num_targets, num_sources, kernel, num_kernel_params, singularity, &
+            approximation, compute_type, degree, max_source_leaf, max_target_leaf, & 
+            verbosity
+
+    DOUBLE PRECISION, POINTER, DIMENSION(:) :: target_x, target_y, target_z, target_q, &
+            source_x, source_y, source_z, source_q, source_w, potential
+
+    DOUBLE PRECISION, DIMENSION(2) :: kernel_params
+
+    DOUBLE PRECISION :: theta, size_check, beta
+
+    CALL MPI_INIT(ierr)
+    CALL MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
+    CALL MPI_COMM_SIZE(MPI_COMM_WORLD, num_proc, ierr)
+
+    kernel = RBS_U
+    num_kernel_params = 1
+    kernel_params(1) = 0.5
+
+    singularity = SKIPPING
+    approximation = LAGRANGE
+    verbosity = 3
+
+    max_source_leaf = 50
+    max_target_leaf = 50
+    size_check = 1.0
+    beta = -1
+
+    compute_type = CLUSTER_CLUSTER
+    theta = 0.7
+    degree = 3
+
+    num_targets = 10000
+    num_sources = 10000
+
+    ALLOCATE(target_x(num_targets), target_y(num_targets), target_z(num_targets), &
+             target_q(num_targets), potential(num_targets), &
+             source_x(num_sources), source_y(num_sources), source_z(num_sources), &
+             source_q(num_sources), source_w(num_sources))
+
+    CALL RANDOM_NUMBER(target_x)
+    CALL RANDOM_NUMBER(target_y)
+    CALL RANDOM_NUMBER(target_z)
+    CALL RANDOM_NUMBER(target_q)
+
+    CALL RANDOM_NUMBER(source_x)
+    CALL RANDOM_NUMBER(source_y)
+    CALL RANDOM_NUMBER(source_z)
+    CALL RANDOM_NUMBER(source_q)
+    CALL RANDOM_NUMBER(source_w)
+
+    ! Calling with kernel as RBS_U
+    CALL BaryTreeInterface(num_targets, num_sources, &
+             C_LOC(target_x), C_LOC(target_y), C_LOC(target_z), C_LOC(target_q), &
+             C_LOC(source_x), C_LOC(source_y), C_LOC(source_z), C_LOC(source_q), &
+             C_LOC(source_w), C_LOC(potential), &
+             kernel, num_kernel_params, C_LOC(kernel_params), &
+             singularity, approximation, compute_type, theta, degree, &
+             max_source_leaf, max_target_leaf, size_check, beta, &
+             verbosity)
+
+    PRINT *, "RBS u total potential is: ", SUM(potential)
+
+    ! Calling with kernel as RBS_V
+    CALL BaryTreeInterface(num_targets, num_sources, &
+             C_LOC(target_x), C_LOC(target_y), C_LOC(target_z), C_LOC(target_q), &
+             C_LOC(source_x), C_LOC(source_y), C_LOC(source_z), C_LOC(source_q), &
+             C_LOC(source_w), C_LOC(potential), &
+             RBS_V, num_kernel_params, C_LOC(kernel_params), &
+             singularity, approximation, compute_type, theta, degree, &
+             max_source_leaf, max_target_leaf, size_check, beta, &
+             verbosity)
+
+    PRINT *, "RBS v total potential is: ", SUM(potential)
+
+    DEALLOCATE(target_x, target_y, target_z, target_q, potential, source_x, source_y, &
+             source_z, source_q, source_w);
+
+    CALL MPI_FINALIZE(ierr)
+
+END PROGRAM

src/CMakeLists.txt

@@ -173,9 +173,26 @@ SET(SRCS_K_MQ
     kernels/mq/mq.h
     kernels/mq/mq_pp.h
     kernels/mq/mq_pp.c
-    kernels/mq/mq.h
+    kernels/mq/mq_pc.h
     kernels/mq/mq_pc.c)
 
+SET(SRCS_K_RBS_U
+    kernels/rbs-u/rbs-u.h
+    kernels/rbs-u/rbs-u_pp.h
+    kernels/rbs-u/rbs-u_pp.c
+    kernels/rbs-u/rbs-u_pc.h
+    kernels/rbs-u/rbs-u_pc.c
+    kernels/rbs-u/rbs-u_cp.h
+    kernels/rbs-u/rbs-u_cp.c)
+
+SET(SRCS_K_RBS_V
+    kernels/rbs-v/rbs-v.h
+    kernels/rbs-v/rbs-v_pp.h
+    kernels/rbs-v/rbs-v_pp.c
+    kernels/rbs-v/rbs-v_pc.h
+    kernels/rbs-v/rbs-v_pc.c
+    kernels/rbs-v/rbs-v_cp.h
+    kernels/rbs-v/rbs-v_cp.c)
 
 SET(SRCS_K_USER
     kernels/user_kernel/user_kernel.h
@@ -195,6 +212,8 @@ SET(SRCS_KERNELS ${SRCS_K_COULOMB}
                  ${SRCS_K_ATAN}
                  ${SRCS_K_SIN_OVER_R}
                  ${SRCS_K_MQ}
+                 ${SRCS_K_RBS_U}
+                 ${SRCS_K_RBS_V}
                  ${SRCS_K_USER})
 
 
@@ -217,6 +236,8 @@ install(TARGETS ${TRGT}                         LIBRARY DESTINATION lib
                                                 ARCHIVE DESTINATION lib)
 install(FILES   interface/BaryTreeInterface.h           DESTINATION include)
 
+install(FILES   ../interfaces/fortran/BaryTreeInterface.fh DESTINATION include)
+
 
 if(ENABLE_GPU_BUILD)
     set(TRGT BaryTree_gpu)
@@ -237,4 +258,6 @@ if(ENABLE_GPU_BUILD)
     install(TARGETS ${TRGT}                       LIBRARY DESTINATION lib
                                                   ARCHIVE DESTINATION lib)
     install(FILES   interface/BaryTreeInterface.h         DESTINATION include)
+
+    install(FILES   ../interfaces/fortran/BaryTreeInterface.fh DESTINATION include)
 endif()