From 5f5fe7dc2895b50714a569d8b644266b12e737e4 Mon Sep 17 00:00:00 2001
From: Abdullah Alperen <alperena@msu.edu>
Date: Fri, 23 Nov 2018 00:23:23 -0500
Subject: [PATCH] PuReMD: fix allocation problem regarding NT buffers

---
 PuReMD/src/allocate.c       |  6 +++---
 PuReMD/src/basic_comm.c     | 32 ++------------------------------
 PuReMD/src/comm_tools.c     |  3 ++-
 PuReMD/src/forces.c         |  6 ++++--
 PuReMD/src/linear_solvers.c | 11 +++--------
 5 files changed, 14 insertions(+), 44 deletions(-)

diff --git a/PuReMD/src/allocate.c b/PuReMD/src/allocate.c
index ffdd965b..d638b8da 100644
--- a/PuReMD/src/allocate.c
+++ b/PuReMD/src/allocate.c
@@ -742,15 +742,15 @@ int  Allocate_MPI_Buffers( mpi_datatypes *mpi_data, int est_recv,
         //TODO: add an exact number
         /* in buffers */
         mpi_data->in_nt_buffer[i] = (void*) 
-            scalloc( my_nbrs[i].est_recv, sizeof(real), "mpibuf:in_nt_buffer", comm );
+            scalloc( my_nt_nbrs[i].est_recv, sizeof(real), "mpibuf:in_nt_buffer", comm );
 
         //TODO
         mpi_buf = &( mpi_data->out_nt_buffers[i] );
         /* allocate storage for the neighbor processor i */
         mpi_buf->index = (int*)
-                         scalloc( my_nbrs[i].est_send, sizeof(int), "mpibuf:nt_index", comm );
+                         scalloc( my_nt_nbrs[i].est_send, sizeof(int), "mpibuf:nt_index", comm );
         mpi_buf->out_atoms = (void*)
-                             scalloc( my_nbrs[i].est_send, sizeof(real), "mpibuf:nt_out_atoms",
+                             scalloc( my_nt_nbrs[i].est_send, sizeof(real), "mpibuf:nt_out_atoms",
                                       comm );
     }
 #endif
diff --git a/PuReMD/src/basic_comm.c b/PuReMD/src/basic_comm.c
index 83f01e52..92c834ef 100644
--- a/PuReMD/src/basic_comm.c
+++ b/PuReMD/src/basic_comm.c
@@ -156,22 +156,11 @@ void Dist_NT( reax_system* system, mpi_datatypes *mpi_data,
             MPI_Irecv( buf + nbr->atoms_str * scale, nbr->atoms_cnt, type,
                     nbr->receive_rank, d, comm, &(req[d]) );
         }
-
-        /*if(d==0)
-        {
-            fprintf( stdout, "p%d dist: entered and need to receive %d numbers from %d\n", system->my_rank, nbr->atoms_cnt, nbr->receive_rank );
-            fflush( stdout );
-        }*/
     }
 
     for( d = 0; d < 6; ++d)
     {
         nbr = &(system->my_nt_nbrs[d]);
-        if(d==0)
-        {
-            //fprintf( stdout, "p%d dist: entered and need to send %d numbers to %d\n", system->my_rank, out_bufs[d].cnt, nbr->rank );
-            //fflush( stdout );
-        }
         /* send both messages in dimension d */
         if ( out_bufs[d].cnt )
         {
@@ -181,16 +170,8 @@ void Dist_NT( reax_system* system, mpi_datatypes *mpi_data,
         }
     }
 
-    /*
-    MPI_Barrier( MPI_COMM_WORLD );
-    fprintf( stdout, "p%d dist: MISSION COMPLETED\n", system->my_rank );
-    fflush( stdout );
-    */
-    
     for( d = 0; d < 6; ++d )
     {
-        //fprintf( stdout, "p%d dist: direction %d\n", system->my_rank, d );
-        //fflush( stdout );
         nbr = &(system->my_nbrs[d]);
         if ( nbr->atoms_cnt )
         {
@@ -198,20 +179,11 @@ void Dist_NT( reax_system* system, mpi_datatypes *mpi_data,
         }
     }
     
-    
     /*for( d = 0; d < count; ++d )
     {
         MPI_Waitany( REAX_MAX_NT_NBRS, req, &index, stat);
     }*/
     
-    
-    /*
-    MPI_Barrier( MPI_COMM_WORLD );
-    fprintf( stdout, "p%d dist: MPI_Waitany is done\n", system->my_rank );
-    fflush( stdout );
-    */
-
-
 #if defined(DEBUG)
     fprintf( stderr, "p%d dist: done\n", system->my_rank );
 #endif
@@ -373,8 +345,8 @@ void Coll_NT( reax_system* system, mpi_datatypes *mpi_data,
     
     for( d = 0; d < 6; ++d )
     {
-        fprintf( stderr, "p%d coll: d = %d\n", system->my_rank, d );
-        fflush( stdout );
+        //fprintf( stderr, "p%d coll: d = %d\n", system->my_rank, d );
+        //fflush( stdout );
         if ( out_bufs[d].cnt )
         {
             MPI_Wait( &(req[d]), &(stat[d]));
diff --git a/PuReMD/src/comm_tools.c b/PuReMD/src/comm_tools.c
index 2162456c..38315557 100644
--- a/PuReMD/src/comm_tools.c
+++ b/PuReMD/src/comm_tools.c
@@ -163,12 +163,13 @@ void Init_Neutral_Territory( reax_system* system, mpi_datatypes *mpi_data )
         if( mpi_data->in_nt_buffer[d] == NULL )
         {
             //TODO
-            mpi_data->in_nt_buffer[d] = (void *) smalloc( 100 * cnt * sizeof(real), "in", comm );
+            mpi_data->in_nt_buffer[d] = (void *) smalloc( SAFER_ZONE * cnt * sizeof(real), "in", comm );
         }
 
         nbr = &(system->my_nt_nbrs[d]);
         nbr->atoms_str = end;
         nbr->atoms_cnt = cnt;
+        nbr->est_recv = MAX( SAFER_ZONE * cnt, MIN_SEND );
         end += cnt;
     }
 }
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 8178856d..93cbe8cf 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -1161,6 +1161,10 @@ void Estimate_Storages( reax_system *system, control_params *control,
                         BO_pi2 = exp( C56 );
                     }
                     else BO_pi2 = C56 = 0.0;
+
+                    /* Initially BO values are the uncorrected ones, page 1 */
+                    BO = BO_s + BO_pi + BO_pi2;
+
                     if ( BO >= control->bo_cut )
                     {
                         ++bond_top[i];
@@ -1278,8 +1282,6 @@ void Compute_Forces( reax_system *system, control_params *control,
     //MPI_Barrier( mpi_data->world );
     if ( system->my_rank == MASTER_NODE )
         Update_Timing_Info( &t_start, &(data->timing.nonb) );
-#endif
-        Update_Timing_Info( &t_start, &(data->timing.nonb) );
 #endif
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "p%d @ step%d: nonbonded forces completed\n",
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index 3606d90d..1d1ccc9a 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -1238,12 +1238,7 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
     else
     {
 #if defined(NEUTRAL_TERRITORY)
-        fprintf(stdout,"BEFORE COLL\n");
-        fflush( stdout );
-
         Coll_NT( system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker );
-        fprintf(stdout,"AFTER COLL\n");
-        fflush( stdout );
 #endif
     }
 
@@ -1288,8 +1283,8 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
 
     for ( i = 0; i < control->cm_solver_max_iters && sqrt(sig_new) / b_norm > tol; ++i )
     {
-        fprintf( stdout, "p%d, i = %d, res =  %.6f\n", system->my_rank, i, sqrt(sig_new) / b_norm );
-        fflush( stdout );
+        //fprintf( stdout, "p%d, i = %d, res =  %.6f\n", system->my_rank, i, sqrt(sig_new) / b_norm );
+        //fflush( stdout );
         t_start = Get_Time( );
 #if defined(NEUTRAL_TERRITORY)
         Dist_NT( system, mpi_data, workspace->d, MPI_DOUBLE, scale, real_packer );
@@ -1364,7 +1359,7 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
         }
 
         t_start = Get_Time( );
-        //sig_new = Parallel_Dot(workspace->r, workspace->p, H->NT, mpi_data->world);
+        sig_old = sig_new;
         sig_new = Parallel_Dot(workspace->r, workspace->p, system->n, mpi_data->world);
         t_allreduce += Get_Timing_Info( t_start );
 
-- 
GitLab