debugging on forcing

include/micm/process/process_set.hpp

@@ -188,20 +188,24 @@ namespace micm
         double rate = cell_rate_constants[i_rxn];
 
         for (std::size_t i_react = 0; i_react < number_of_reactants_[i_rxn]; ++i_react)
-          rate *= cell_state[react_id[i_react]];
-        for (std::size_t i_react = 0; i_react < number_of_reactants_[i_rxn]; ++i_react){
-          std::cout << "this is cell forcing before update: " << cell_forcing[react_id[i_react]] <<std::endl; 
+          rate *= cell_state[react_id[i_react]]; 
+
+        for (std::size_t i_react = 0; i_react < number_of_reactants_[i_rxn]; ++i_react)
           cell_forcing[react_id[i_react]] -= rate;
-          std::cout << "this is cell forcing after update: "  << cell_forcing[react_id[i_react]] <<std::endl; 
-          }
 
-        for (std::size_t i_prod = 0; i_prod < number_of_products_[i_rxn]; ++i_prod)
-          cell_forcing[prod_id[i_prod]] += yield[i_prod] * rate;
+
+        // for (std::size_t i_prod = 0; i_prod < number_of_products_[i_rxn]; ++i_prod)
+        //   cell_forcing[prod_id[i_prod]] += yield[i_prod] * rate;
+
         react_id += number_of_reactants_[i_rxn];
         prod_id += number_of_products_[i_rxn];
         yield += number_of_products_[i_rxn];
       }
     }
+    double* forcing_data = forcing.AsVector().data(); 
+        for (int i = 0; i < forcing.AsVector().size(); i++){
+          std::cout << "this is cell forcing after update: "<< forcing_data[i]<<std::endl; 
+        }
   };
 
   inline void ProcessSet::AddJacobianTerms(

src/process/process_set.cu

@@ -55,9 +55,9 @@ namespace micm {
             int state_forcing_col_index = reactant_ids_[reactant_ids_index]; 
 
             forcing[row_index * state_forcing_columns + state_forcing_col_index] -=rate; 
-            if (tid == 1){
-                cell_forcing[tid + i_reactant] = forcing[row_index * state_forcing_columns + state_forcing_col_index]; 
-            }
+
+                cell_forcing[row_index * state_forcing_columns + state_forcing_col_index] = forcing[row_index * state_forcing_columns + state_forcing_col_index]; 
+
         }
 
         for (int i_product = 0; i_product < product_num; i_product++){
@@ -143,9 +143,10 @@ namespace micm {
 
         double* d_cell_forcing; 
         double* cell_forcing; 
-        cudaMalloc(&d_cell_forcing, sizeof(double) * 8); 
-        cell_forcing = (double*)malloc(sizeof(double) * 8); 
-
+        cudaMalloc(&d_cell_forcing, sizeof(double) * 10); 
+        cell_forcing = (double*)malloc(sizeof(double) * 10); 
+        cudaMemcpy(d_cell_forcing, forcing_data, sizeof(double)* 10, cudaMemcpyHostToDevice); 
+
 
 
         //allocate device memory
@@ -209,26 +210,26 @@ namespace micm {
 
         cudaMemcpy(forcing_data, d_forcing, state_forcing_bytes, cudaMemcpyDeviceToHost);
 
-        // //debugging 
-        // cudaMemcpy(rate_array, d_rate_array, sizeof(double)*rate_array_size, cudaMemcpyDeviceToHost);    
-        // std::cout << "this is rate_array before update: "<< std::endl; 
-        // for (int k = 0; k < rate_array_size; k++){
-        //     std::cout << rate_array[k]<<std::endl; 
-        // }
-        // cudaMemcpy(state_variable, d_state_variable, sizeof(double)*2, cudaMemcpyDeviceToHost);   
-        // std::cout << "This is state variable for first thread"<<std::endl;
-        // for (int k = 0; k < 2; k++){
-        //     std::cout << state_variable[k]<<std::endl; 
-        // }
+        //debugging 
+        cudaMemcpy(rate_array, d_rate_array, sizeof(double)*rate_array_size, cudaMemcpyDeviceToHost);    
+        std::cout << "this is rate_array before update: "<< std::endl; 
+        for (int k = 0; k < rate_array_size; k++){
+            std::cout << rate_array[k]<<std::endl; 
+        }
+        cudaMemcpy(state_variable, d_state_variable, sizeof(double)*2, cudaMemcpyDeviceToHost);   
+        std::cout << "This is state variable for first thread"<<std::endl;
+        for (int k = 0; k < 2; k++){
+            std::cout << state_variable[k]<<std::endl; 
+        }
 
-        // cudaMemcpy(rate_array_post, d_rate_array_post, sizeof(double)*rate_array_size, cudaMemcpyDeviceToHost);    
-        // std::cout << "this is rate_array after update: "<< std::endl; 
-        // for (int k = 0; k < rate_array_size; k++){
-        //     std::cout << rate_array_post[k]<<std::endl; 
-        // }
+        cudaMemcpy(rate_array_post, d_rate_array_post, sizeof(double)*rate_array_size, cudaMemcpyDeviceToHost);    
+        std::cout << "this is rate_array after update: "<< std::endl; 
+        for (int k = 0; k < rate_array_size; k++){
+            std::cout << rate_array_post[k]<<std::endl; 
+        }
 
-       cudaMemcpy(cell_forcing, d_cell_forcing, sizeof(double) * 2, cudaMemcpyDeviceToHost); 
-       for (int k = 0; k < 8; k++){
+       cudaMemcpy(cell_forcing, d_cell_forcing, sizeof(double)*10, cudaMemcpyDeviceToHost); 
+       for (int k = 0; k < 10; k++){
         std::cout << "this is cell forcing after update"<< cell_forcing[k]<<std::endl; 
        }
 
@@ -253,4 +254,4 @@ namespace micm {
 
         }
     }//namespace cuda 
-}//namespace micm
+}//namespace micm

test/unit/process/test_cuda_process_set.cpp

@@ -59,6 +59,10 @@ TEST(ProcessSet, Constructor)
   int product_ids_size = set.product_ids_vector().size();
   const double* yields = set.yields_vector().data();
   int yields_size = set.yields_vector().size(); 
+
+  for (int i = 0; i < number_of_reactants_size; i++){
+    std::cout << number_of_reactants[i]<<std::endl; 
+  }
 
   micm::cuda::AddForcingTerms_kernelSetup(
     number_of_reactants,

test/unit/process/test_process_set.cpp

@@ -47,7 +47,7 @@ TEST(ProcessSet, Constructor)
   rate_constants[1] = { 110.0, 120.0, 130.0 };
 
   micm::Matrix<double> forcing{ 2, 5, 1000.0 };
-
+  
   set.AddForcingTerms(rate_constants, state.variables_, forcing);
 
   EXPECT_EQ(forcing[0][0], 1000.0 - 10.0 * 0.1 * 0.3 + 20.0 * 0.2);