diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index 8420f82..2cf6144 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -70,6 +70,7 @@ sink_gravity(int3 globalVertexIdx){
 
 
 #if LSINK
+// Give Truelove density
 Scalar
 truelove_density(in ScalarField lnrho){
     const Scalar rho = exp(value(lnrho));
@@ -77,10 +78,7 @@ truelove_density(in ScalarField lnrho){
     const Scalar TJ_rho = ((M_PI) * ((dsx * dsx) / Jeans_length_squared) * cs2_sound) / (AC_G_const * dsx * dsx);
     //TODO: dsx will cancel out, deal with it later for optimization.      
  
-    Scalar accretion_rho = rho - TJ_rho;
-    if (accretion_rho < 0){
-        accretion_rho = Scalar(0);
-    } 
+    Scalar accretion_rho = TJ_rho;
        
     return accretion_rho;
 }
@@ -96,16 +94,7 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
     const Scalar profile_range = DCONST_REAL(AC_accretion_range);
     const Scalar accretion_distance = length(grid_pos - sink_pos);   
     Scalar accretion_density; 
-//    if ((accretion_distance) <= profile_range){
-//        // calculate accretion according to chosen criterion for the grid cell.
-//        accretion_density = truelove_density(lnrho);
-//    } else {
-//        accretion_density = Scalar(0.0); 
-//    }
-//    return accretion_density;
      
-    // multiplying the truelove density by a wave function to avoid step-function like accretion profile.   
-
     Scalar weight;
 //    const Scalar weight = exp(-(accretion_distance/profile_range));   
     // Step function weighting 
@@ -115,17 +104,16 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
         weight = Scalar(0.0);
     }
 
-    const Scalar lnrho_min = Scalar(-10.0);  //TODO Define from astaroth.conf
+//    const Scalar lnrho_min = Scalar(-10.0);  //TODO Define from astaroth.conf
+    const Scalar lnrho_min = log(truelove_density(lnrho));
 //    const Scalar sink_mass = DCONST_REAL(AC_M_sink);
-//    const Scalar rate = truelove_density(lnrho);
 //    const Scalar B = Scalar(0.5);
 //    const Scalar k = Scalar(1.5);
 //    const Scalar rate = B * (pow(sink_mass, k) / (dsx * dsy * dsz));
     Scalar rate;     
     if (value(lnrho) > lnrho_min) {
-        rate = exp(value(lnrho)) / dt;
+        rate = (exp(value(lnrho)) - exp(lnrho_min)) / dt;
     } else { 
-//        rate = exp(lnrho_min) / dt;
         rate = Scalar(0.0);
     }
     accretion_density = weight * rate ;
diff --git a/config/templates/sinktest.conf b/config/templates/sinktest.conf
index 63ac4b8..258e0a3 100644
--- a/config/templates/sinktest.conf
+++ b/config/templates/sinktest.conf
@@ -54,7 +54,7 @@ AC_sink_pos_x = 3.14
 AC_sink_pos_y = 3.14
 AC_sink_pos_z = 3.14
 //AC_M_sink_Msun = 0.005
-AC_M_sink_Msun = 1e-6
+AC_M_sink_Msun = 1e-8
 AC_soft = 0.36
 
 // Accretion Parameters
diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc
index cf5d05e..c883086 100644
--- a/src/standalone/simulation.cc
+++ b/src/standalone/simulation.cc
@@ -250,9 +250,13 @@ run_simulation(void)
         const AcReal dt   = host_timestep(umax, mesh_info);
 
 #if LSINK
+
         const AcReal sum_mass = acReduceScal(RTYPE_MAX, VTXBUF_ACCRETION);
 	accreted_mass = accreted_mass + sum_mass;
-        AcReal sink_mass = AC_M_sink_init + accreted_mass;
+        AcReal sink_mass = 0.0;
+        //if (i > 1000 ) {
+	    sink_mass = AC_M_sink_init + accreted_mass;
+	//} 
         printf("sink mass is: %e \n", sink_mass); 
         printf("accreted mass is: %e \n", accreted_mass); 
         acLoadDeviceConstant(AC_M_sink, sink_mass);