diff --git a/acc/mhd_solver/stencil_definition.sdh b/acc/mhd_solver/stencil_definition.sdh
index c4324e5..abc341c 100644
--- a/acc/mhd_solver/stencil_definition.sdh
+++ b/acc/mhd_solver/stencil_definition.sdh
@@ -5,6 +5,7 @@
 #define LTEMPERATURE (0)
 #define LFORCING (1)
 #define LUPWD (1)
+#define LSINK (0)
 
 #define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter
 
@@ -37,6 +38,16 @@ uniform Scalar AC_star_pos_x;
 uniform Scalar AC_star_pos_y;
 uniform Scalar AC_star_pos_z;
 uniform Scalar AC_M_star;
+// properties of sink particle
+uniform Scalar AC_sink_pos_x;
+uniform Scalar AC_sink_pos_y;
+uniform Scalar AC_sink_pos_z;
+uniform Scalar AC_M_sink;
+uniform Scalar AC_M_sink_init;
+uniform Scalar AC_M_sink_Msun;
+uniform Scalar AC_soft;
+uniform Scalar AC_accretion_range;
+uniform Scalar AC_switch_accretion;
 //  Run params
 uniform Scalar AC_cdt;
 uniform Scalar AC_cdtv;
@@ -78,8 +89,9 @@ uniform Scalar AC_ff_hel_imx;
 uniform Scalar AC_ff_hel_imy;
 uniform Scalar AC_ff_hel_imz;
 //  Additional helper params  //  (deduced from other params do not set these directly!)
-uniform Scalar AC_G_CONST;
+uniform Scalar AC_G_const;
 uniform Scalar AC_GM_star;
+uniform Scalar AC_unit_mass;
 uniform Scalar AC_sq2GM_star;
 uniform Scalar AC_cs2_sound;
 uniform Scalar AC_inv_dsx;
@@ -116,3 +128,8 @@ uniform ScalarField VTXBUF_UUZ;
 #else
 uniform ScalarField VTXBUF_LNRHO;
 #endif
+
+#if LSINK
+uniform ScalarField VTXBUF_ACCRETION;
+#endif
+
diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index 0db62e0..a0787c7 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -23,20 +23,155 @@ gradients(in VectorField uu)
     return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)};
 }
 
+#if LSINK 
+Vector
+sink_gravity(int3 globalVertexIdx){
+    int accretion_switch = int(AC_switch_accretion);
+    if (accretion_switch == 1){
+        Vector force_gravity;
+        const Vector grid_pos = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
+                                         (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
+                                         (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz};
+        const Scalar sink_mass = AC_M_sink;
+        const Vector sink_pos = (Vector){AC_sink_pos_x,
+                                         AC_sink_pos_y,
+                                         AC_sink_pos_z}; 
+        const Scalar distance = length(grid_pos - sink_pos);
+        const Scalar soft = AC_soft;
+        //MV: The commit 083ff59 had AC_G_const defined wrong here in DSL making it exxessively strong.
+        //MV: Scalar gravity_magnitude = ... below is correct! 
+        const Scalar gravity_magnitude = (AC_G_const * sink_mass) / pow(((distance * distance) +  soft*soft), 1.5);
+        const Vector direction = (Vector){(sink_pos.x - grid_pos.x) / distance,
+                                          (sink_pos.y - grid_pos.y) / distance,
+                                          (sink_pos.z - grid_pos.z) / distance};
+        force_gravity = gravity_magnitude * direction;
+        return force_gravity;
+    } else {
+        return (Vector){0.0, 0.0, 0.0};
+    }
+}
+#endif
+
+
+#if LSINK
+// Give Truelove density
 Scalar
-continuity(in VectorField uu, in ScalarField lnrho)
+truelove_density(in ScalarField lnrho){
+    const Scalar rho = exp(value(lnrho));
+    const Scalar Jeans_length_squared = (M_PI * AC_cs2_sound) / (AC_G_const * rho);
+    const Scalar TJ_rho = ((M_PI) * ((AC_dsx * AC_dsx) / Jeans_length_squared) * AC_cs2_sound) / (AC_G_const * AC_dsx * AC_dsx);
+    //TODO: AC_dsx will cancel out, deal with it later for optimization.      
+ 
+    Scalar accretion_rho = TJ_rho;
+       
+    return accretion_rho;
+}
+
+// This controls accretion of density/mass to the sink particle. 
+Scalar
+sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
+    const Vector grid_pos = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
+                                     (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
+                                     (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz};
+    const Vector sink_pos = (Vector){AC_sink_pos_x,
+                                     AC_sink_pos_y,
+                                     AC_sink_pos_z};
+    const Scalar profile_range = AC_accretion_range;
+    const Scalar accretion_distance = length(grid_pos - sink_pos);
+    int accretion_switch = AC_switch_accretion; 
+    Scalar accretion_density; 
+    Scalar weight;
+
+    if (accretion_switch == 1){   
+        if ((accretion_distance) <= profile_range){
+            //weight = Scalar(1.0);
+            //Hann window function
+            Scalar window_ratio = accretion_distance/profile_range;
+            weight = Scalar(0.5)*(Scalar(1.0) - cos(Scalar(2.0)*M_PI*window_ratio));
+        } else {
+            weight = Scalar(0.0);
+        }
+        
+        //Truelove criterion is used as a kind of arbitrary density floor. 
+        const Scalar lnrho_min = log(truelove_density(lnrho));
+        Scalar rate;     
+        if (value(lnrho) > lnrho_min) {
+            rate = (exp(value(lnrho)) - exp(lnrho_min)) / dt;
+        } else { 
+            rate = Scalar(0.0);
+        }
+        accretion_density = weight * rate ;
+    } else {
+        accretion_density = Scalar(0.0); 
+    }
+    return accretion_density;
+} 
+
+// This controls accretion of velocity to the sink particle. 
+Vector
+sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) {
+    const Vector grid_pos = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
+                                     (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
+                                     (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz};
+    const Vector sink_pos = (Vector){AC_sink_pos_x,
+                                     AC_sink_pos_y,
+                                     AC_sink_pos_z};
+    const Scalar profile_range = AC_accretion_range;
+    const Scalar accretion_distance = length(grid_pos - sink_pos);   
+    int accretion_switch = AC_switch_accretion; 
+    Vector accretion_velocity; 
+
+    if (accretion_switch == 1){
+        Scalar weight;
+        // Step function weighting
+        // Arch of a cosine function? 
+        // Cubic spline x^3 - x in range [-0.5 , 0.5] 
+        if ((accretion_distance) <= profile_range){
+            //weight = Scalar(1.0);
+            //Hann window function
+            Scalar window_ratio = accretion_distance/profile_range;
+            weight = Scalar(0.5)*(Scalar(1.0) - cos(Scalar(2.0)*M_PI*window_ratio));
+        } else {
+            weight = Scalar(0.0);
+        }
+
+
+        Vector rate;     
+        // MV: Could we use divergence here ephasize velocitie which are compressive and 
+        // MV: not absorbins stuff that would not be accreted anyway? 
+        if (length(value(uu)) > Scalar(0.0)) {
+            rate = (Scalar(1.0)/dt) * value(uu);
+        } else { 
+            rate = (Vector){0.0, 0.0, 0.0};
+        }
+        accretion_velocity = weight * rate ;
+    } else {
+        accretion_velocity = (Vector){0.0, 0.0, 0.0};
+    }
+    return accretion_velocity;
+}
+#endif
+
+
+Scalar
+continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) 
 {
     return -dot(value(uu), gradient(lnrho))
 #if LUPWD
            // This is a corrective hyperdiffusion term for upwinding.
            + upwd_der6(uu, lnrho)
+#endif
+#if LSINK
+	   - sink_accretion(globalVertexIdx, lnrho, dt) / exp(value(lnrho)) 
 #endif
            - divergence(uu);
 }
 
+
+
 #if LENTROPY
 Vector
-momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa)
+momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa, Scalar dt) 
 {
     const Matrix S   = stress_tensor(uu);
     const Scalar cs2 = AC_cs2_sound * exp(AC_gamma * value(ss) / AC_cp_sound +
@@ -55,12 +190,22 @@ momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorFi
                        AC_nu_visc *
                            (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
                             Scalar(2.) * mul(S, gradient(lnrho))) +
-                       AC_zeta * gradient_of_divergence(uu);
+                       AC_zeta * gradient_of_divergence(uu)
+    #if LSINK
+                       //Gravity term
+                       + sink_gravity(globalVertexIdx) 
+                       //Corresponding loss of momentum
+                       -     //(Scalar(1.0) / Scalar( (AC_dsx*AC_dsy*AC_dsz) * exp(value(lnrho)))) *  // Correction factor by unit mass
+    	                 sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014)
+                       ; 
+    #else
+                       ;
+    #endif
     return mom;
 }
 #elif LTEMPERATURE
 Vector
-momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
+momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in ScalarField tt) 
 {
     Vector mom;
 
@@ -72,17 +217,21 @@ momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
     mom = -mul(gradients(uu), value(uu)) - pressure_term +
           AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
                         Scalar(2.) * mul(S, gradient(lnrho))) +
-          AC_zeta * gradient_of_divergence(uu);
+          AC_zeta * gradient_of_divergence(uu)
+    #if LSINK
+          + sink_gravity(globalVertexIdx);
+    #else
+                                         ;
+    #endif
 
 #if LGRAVITY
     mom = mom - (Vector){0, 0, -10.0};
 #endif
-
     return mom;
 }
 #else
 Vector
-momentum(in VectorField uu, in ScalarField lnrho)
+momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) 
 {
     Vector mom;
 
@@ -93,7 +242,16 @@ momentum(in VectorField uu, in ScalarField lnrho)
     mom = -mul(gradients(uu), value(uu)) - AC_cs2_sound * gradient(lnrho) +
           AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
                         Scalar(2.) * mul(S, gradient(lnrho))) +
-          AC_zeta * gradient_of_divergence(uu);
+          AC_zeta * gradient_of_divergence(uu)
+    #if LSINK
+          + sink_gravity(globalVertexIdx)
+          //Corresponding loss of momentum
+          -     //(Scalar(1.0) / Scalar( (AC_dsx*AC_dsy*AC_dsz) * exp(value(lnrho)))) *  // Correction factor by unit mass
+    	    sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014)
+                                              ; 
+    #else
+                                              ;
+    #endif
 
 #if LGRAVITY
     mom = mom - (Vector){0, 0, -10.0};
@@ -184,15 +342,23 @@ heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
 
 #if LFORCING
 Vector
-simple_vortex_forcing(Vector a, Vector b, Scalar magnitude)
-{
-    return magnitude * cross(normalized(b - a), (Vector){0, 0, 1}); // Vortex
-}
-
+    simple_vortex_forcing(Vector a, Vector b, Scalar magnitude){
+    int accretion_switch = AC_switch_accretion; 
+    
+    if (accretion_switch == 0){
+        return magnitude * cross(normalized(b - a), (Vector){ 0, 0, 1}); // Vortex
+    } else { 
+        return (Vector){0,0,0};  
+    }
+}        
 Vector
-simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude)
-{
-    return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
+    simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude){
+    int accretion_switch = AC_switch_accretion;
+    if (accretion_switch == 0){            
+        return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
+    } else {
+        return (Vector){0,0,0};
+    }
 }
 
 // The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable
@@ -234,38 +400,42 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto
 Vector
 forcing(int3 globalVertexIdx, Scalar dt)
 {
-    Vector a         = Scalar(.5) * (Vector){globalGridN.x * AC_dsx, globalGridN.y * AC_dsy,
-                                     globalGridN.z * AC_dsz}; // source (origin)
-    Vector xx        = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
-                         (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
-                         (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; // sink (current index)
-    const Scalar cs2 = AC_cs2_sound;
-    const Scalar cs  = sqrt(cs2);
+    int accretion_switch = AC_switch_accretion;
+    if (accretion_switch == 0){
 
-    // Placeholders until determined properly
-    Scalar magnitude = AC_forcing_magnitude;
-    Scalar phase     = AC_forcing_phase;
-    Vector k_force   = (Vector){AC_k_forcex, AC_k_forcey, AC_k_forcez};
-    Vector ff_re     = (Vector){AC_ff_hel_rex, AC_ff_hel_rey, AC_ff_hel_rez};
-    Vector ff_im     = (Vector){AC_ff_hel_imx, AC_ff_hel_imy, AC_ff_hel_imz};
-
-    // Determine that forcing funtion type at this point.
-    // Vector force = simple_vortex_forcing(a, xx, magnitude);
-    // Vector force = simple_outward_flow_forcing(a, xx, magnitude);
-    Vector force = helical_forcing(magnitude, k_force, xx, ff_re, ff_im, phase);
-
-    // Scaling N = magnitude*cs*sqrt(k*cs/dt)  * dt
-    const Scalar NN = cs * sqrt(AC_kaver * cs);
-    // MV: Like in the Pencil Code. I don't understandf the logic here.
-    force.x = sqrt(dt) * NN * force.x;
-    force.y = sqrt(dt) * NN * force.y;
-    force.z = sqrt(dt) * NN * force.z;
-
-    if (is_valid(force)) {
-        return force;
-    }
-    else {
-        return (Vector){0, 0, 0};
+        Vector a = Scalar(.5) * (Vector){globalGridN.x * AC_dsx,
+                                         globalGridN.y * AC_dsy,
+                                         globalGridN.z * AC_dsz}; // source (origin)
+        Vector xx = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
+                             (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
+                             (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; // sink (current index)
+        const Scalar cs2 = AC_cs2_sound;
+        const Scalar cs = sqrt(cs2);
+    
+        //Placeholders until determined properly
+        Scalar magnitude = AC_forcing_magnitude;
+        Scalar phase     = AC_forcing_phase;
+        Vector k_force   = (Vector){AC_k_forcex,   AC_k_forcey,   AC_k_forcez};
+        Vector ff_re     = (Vector){AC_ff_hel_rex, AC_ff_hel_rey, AC_ff_hel_rez};
+        Vector ff_im     = (Vector){AC_ff_hel_imx, AC_ff_hel_imy, AC_ff_hel_imz};
+    
+    
+        //Determine that forcing funtion type at this point.
+        //Vector force = simple_vortex_forcing(a, xx, magnitude);
+        //Vector force = simple_outward_flow_forcing(a, xx, magnitude);
+        Vector force   = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase);
+    
+        //Scaling N = magnitude*cs*sqrt(k*cs/dt)  * dt
+        const Scalar NN = cs*sqrt(AC_kaver*cs);
+        //MV: Like in the Pencil Code. I don't understandf the logic here.
+        force.x = sqrt(dt)*NN*force.x;
+        force.y = sqrt(dt)*NN*force.y;
+        force.z = sqrt(dt)*NN*force.z;
+    
+        if (is_valid(force)) { return force; }
+        else                 { return (Vector){0, 0, 0}; }
+    } else {
+          return (Vector){0,0,0};
     }
 }
 #endif // LFORCING
@@ -293,24 +463,29 @@ in ScalarField tt(VTXBUF_TEMPERATURE);
 out ScalarField out_tt(VTXBUF_TEMPERATURE);
 #endif
 
+#if LSINK
+in ScalarField accretion(VTXBUF_ACCRETION);
+out ScalarField out_accretion(VTXBUF_ACCRETION);
+#endif
+
 Kernel void
 solve()
 {
     Scalar dt = AC_dt;
-    out_lnrho = rk3(out_lnrho, lnrho, continuity(uu, lnrho), dt);
+    out_lnrho = rk3(out_lnrho, lnrho, continuity(globalVertexIdx, uu, lnrho, dt), dt);
 
 #if LMAGNETIC
     out_aa = rk3(out_aa, aa, induction(uu, aa), dt);
 #endif
 
 #if LENTROPY
-    out_uu = rk3(out_uu, uu, momentum(uu, lnrho, ss, aa), dt);
+    out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho, ss, aa, dt), dt);
     out_ss = rk3(out_ss, ss, entropy(ss, uu, lnrho, aa), dt);
 #elif LTEMPERATURE
-    out_uu = rk3(out_uu, uu, momentum(uu, lnrho, tt), dt);
+    out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho, tt, dt), dt);
     out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt);
 #else
-    out_uu = rk3(out_uu, uu, momentum(uu, lnrho), dt);
+    out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho, dt), dt);
 #endif
 
 #if LFORCING
@@ -318,4 +493,12 @@ solve()
         out_uu = out_uu + forcing(globalVertexIdx, dt);
     }
 #endif
+
+#if LSINK
+    out_accretion = rk3(out_accretion, accretion, sink_accretion(globalVertexIdx, lnrho, dt), dt);// unit now is rho!
+          
+    if (step_number == 2) {
+        out_accretion = out_accretion * AC_dsx * AC_dsy * AC_dsz;// unit is now mass!
+    }
+#endif
 }
diff --git a/analysis/python/astar/data/read.py b/analysis/python/astar/data/read.py
index 270263f..f55f145 100644
--- a/analysis/python/astar/data/read.py
+++ b/analysis/python/astar/data/read.py
@@ -78,6 +78,8 @@ class Mesh:
         if self.ok:
 
             self.ss, timestamp, ok = read_bin('VTXBUF_ENTROPY', fdir, fnum, self.minfo)
+
+            self.accretion, timestamp, ok = read_bin('VTXBUF_ACCRETION', fdir, fnum, self.minfo)
  
             #TODO Generalize is a dict. Do not hardcode!  
             uux, timestamp, ok = read_bin('VTXBUF_UUX', fdir, fnum, self.minfo)
@@ -116,6 +118,7 @@ def parse_ts(fdir, fname):
         line[i] = line[i].replace('VTXBUF_', "")
         line[i] = line[i].replace('UU', "uu")
         line[i] = line[i].replace('_total', "tot")
+        line[i] = line[i].replace('ACCRETION', "acc")
         line[i] = line[i].replace('A', "aa")
         line[i] = line[i].replace('LNRHO', "lnrho")
         line[i] = line[i].replace('ENTROPY', "ss")
diff --git a/analysis/python/astar/visual/lineplot.py b/analysis/python/astar/visual/lineplot.py
index f33d75d..45d39df 100644
--- a/analysis/python/astar/visual/lineplot.py
+++ b/analysis/python/astar/visual/lineplot.py
@@ -32,7 +32,7 @@ def plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3):
     plt.xlabel(xaxis)
     plt.legend()
 
-def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False, uuz=False, ss=False):
+def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False, uuz=False, ss=False, acc=False, sink=False):
 
     if show_all:
         lnrho=True
@@ -41,6 +41,8 @@ def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False,
         uuy=True
         uuz=True
         ss=True
+        acc=True
+        sink=True
 
     if lnrho:
         plt.figure()
@@ -89,5 +91,27 @@ def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False,
         yaxis2 = 'ss_min'
         yaxis3 = 'ss_max'
         plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
+
+    if acc:   
+        plt.figure()
+        xaxis = 't_step'
+        yaxis1 = 'acc_rms'
+        yaxis2 = 'acc_min'
+        yaxis3 = 'acc_max'
+        plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
+
+        plt.figure()
+        xaxis = 't_step'
+        yaxis1 = 'sink_mass'
+        yaxis2 = 'sink_mass'
+        yaxis3 = 'sink_mass'
+        plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
+  
+        plt.figure()
+        xaxis = 't_step'
+        yaxis1 = 'accreted_mass'
+        yaxis2 = 'accreted_mass'
+        yaxis3 = 'accreted_mass'
+        plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
   
     plt.show()
diff --git a/analysis/python/jupyter/notebook_example.ipynb b/analysis/python/jupyter/notebook_example.ipynb
index f7e4956..3cbc3ce 100644
--- a/analysis/python/jupyter/notebook_example.ipynb
+++ b/analysis/python/jupyter/notebook_example.ipynb
@@ -24,8 +24,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "meshdir  = \"/scratch/data/mvaisala/forcingtest/\"\n",
-    "#meshdir  = \"/scratch/data/mvaisala/normaltest/\""
+    "meshdir  = \"/home/mvaisala/astaroth/build/\""
    ]
   },
   {
@@ -35,7 +34,7 @@
    "outputs": [],
    "source": [
     "#imesh = 30000\n",
-    "imesh = 100\n",
+    "imesh = 0\n",
     "mesh = ad.read.Mesh(imesh, fdir=meshdir)"
    ]
   },
@@ -74,10 +73,12 @@
     "vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$')\n",
     "vis.slices.plot_3(mesh, mesh.uu[1], title = r'$u_y$')\n",
     "vis.slices.plot_3(mesh, mesh.uu[2], title = r'$u_z$')\n",
-    "vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\\ln_\\rho$')\n",
+    "vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\\ln_\\rho$')#, colrange=[0,0.1])\n",
     "vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$\\rho$')\n",
-    "vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$N_\\mathrm{col}$', slicetype = 'sum')\n",
-    "vis.slices.plot_3(mesh, mesh.ss, title = r'$s$')\n"
+    "vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$N_\\mathrm{col}$', slicetype = 'sum')#, colrange=[130,139])\n",
+    "vis.slices.plot_3(mesh, uu_tot, title = r'$\\Sigma \\|u\\|$', slicetype = 'sum')#, colrange=[130,139])\n",
+    "#vis.slices.plot_3(mesh, mesh.ss, title = r'$s$')\n",
+    "vis.slices.plot_3(mesh, mesh.accretion, title = r'$Accretion buffer$')\n"
    ]
   },
   {
@@ -102,8 +103,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, uuy=1, uuz=1, ss=1)"
+    "vis.lineplot.plot_ts(ts, lnrho=True, acc=True, uux=True)\n",
+    "\n"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -122,7 +131,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.2"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,
diff --git a/analysis/python/samples/readtest.py b/analysis/python/samples/readtest.py
index b2e4458..7f22858 100644
--- a/analysis/python/samples/readtest.py
+++ b/analysis/python/samples/readtest.py
@@ -45,7 +45,7 @@ print("sys.argv", sys.argv)
 #meshdir = "/tiara/home/mvaisala/astaroth-code/astaroth_2.0/build/"
 #meshdir  = "/tiara/ara/data/mvaisala/tmp/astaroth-code/astaroth_2.0/build/"
 #meshdir = "/tiara/ara/data/mvaisala/asth_testbed_double/"
-meshdir  = "/scratch/data/mvaisala/forcingtest/"
+meshdir  = "/home/mvaisala/astaroth/build/"
 
 if "xtopbound" in sys.argv: 
     for i in range(0, 171):
@@ -199,6 +199,6 @@ if 'sl' in sys.argv:
 
 if 'ts' in sys.argv:
    ts = ad.read.TimeSeries(fdir=meshdir)
-   vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, ss=1)
+   vis.lineplot.plot_ts(ts, show_all=True)
 
 
diff --git a/config/astaroth.conf b/config/astaroth.conf
index 32f50a3..a003c8b 100644
--- a/config/astaroth.conf
+++ b/config/astaroth.conf
@@ -49,6 +49,22 @@ AC_gamma = 0.5
 AC_lnT0 = 1.2
 AC_lnrho0 = 1.3
 
+// Sink Particle
+AC_sink_pos_x = 3.14
+AC_sink_pos_y = 3.14
+AC_sink_pos_z = 3.14
+AC_M_sink_Msun = 1.0
+AC_soft = 0.12
+
+// Accretion Parameters
+// profile_range is multiple of dsx
+AC_accretion_range = 2.0
+
+// Physical properties of the domain
+AC_unit_velocity = 1.0
+AC_unit_density = 1.0
+AC_unit_length = 1.0
+
 /*
  * =============================================================================
  * Initial conditions
diff --git a/config/templates/core_collapse.conf b/config/templates/core_collapse.conf
new file mode 100644
index 0000000..df70dbe
--- /dev/null
+++ b/config/templates/core_collapse.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 128
+AC_ny = 128
+AC_nz = 128
+
+AC_dsx = 0.04908738521
+AC_dsy = 0.04908738521
+AC_dsz = 0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 20001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 2e10
diff --git a/config/templates/sink_aftersummer/128_basic_rotating.conf b/config/templates/sink_aftersummer/128_basic_rotating.conf
new file mode 100644
index 0000000..ccfc07c
--- /dev/null
+++ b/config/templates/sink_aftersummer/128_basic_rotating.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 128
+AC_ny = 128
+AC_nz = 128
+
+AC_dsx = 0.04908738521
+AC_dsy = 0.04908738521
+AC_dsz = 0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 501
+AC_save_steps = 1
+AC_bin_steps = 100
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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 = 1.0 // 1e-4 //1e-5 
+AC_soft = 0.25
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 1.2
+AC_ampl_uu = 1.0
diff --git a/config/templates/sinktest.conf b/config/templates/sinktest.conf
new file mode 100644
index 0000000..b1e118e
--- /dev/null
+++ b/config/templates/sinktest.conf
@@ -0,0 +1,77 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 128
+AC_ny = 128
+AC_nz = 128
+
+AC_dsx = 0.04908738521
+AC_dsy = 0.04908738521
+AC_dsz = 0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 1001
+AC_save_steps = 1
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-3
+AC_cs_sound = 1.0
+AC_zeta = 0.01
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-5
+AC_kmin              = 0.8
+AC_kmax              = 1.2 
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-10
+AC_soft = 0.36
+
+// Accretion Parameters
+// profile_range shoul be close to a multiple of dsx. E.g. 4*dsx
+AC_accretion_range = 0.2
+
+// Physical properties of the domain
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3
+AC_unit_density = 3e-22
+// using 10,000 AU 
+AC_unit_length = 1.5e17
+
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 0.0
diff --git a/config/templates/sinktest/128res_non-rotating_highvisc.conf b/config/templates/sinktest/128res_non-rotating_highvisc.conf
new file mode 100644
index 0000000..3a798dc
--- /dev/null
+++ b/config/templates/sinktest/128res_non-rotating_highvisc.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 128
+AC_ny = 128
+AC_nz = 128
+
+AC_dsx = 0.04908738521
+AC_dsy = 0.04908738521
+AC_dsz = 0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 20001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 0.0
diff --git a/config/templates/sinktest/128res_rotating_highvisc.conf b/config/templates/sinktest/128res_rotating_highvisc.conf
new file mode 100644
index 0000000..df70dbe
--- /dev/null
+++ b/config/templates/sinktest/128res_rotating_highvisc.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 128
+AC_ny = 128
+AC_nz = 128
+
+AC_dsx = 0.04908738521
+AC_dsy = 0.04908738521
+AC_dsz = 0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 20001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 2e10
diff --git a/config/templates/sinktest/256res_non-rotating_highvisc.conf b/config/templates/sinktest/256res_non-rotating_highvisc.conf
new file mode 100644
index 0000000..e451e7b
--- /dev/null
+++ b/config/templates/sinktest/256res_non-rotating_highvisc.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 256
+AC_ny = 256 
+AC_nz = 256
+
+AC_dsx = 0.0245436926 //0.04908738521
+AC_dsy = 0.0245436926 //0.04908738521
+AC_dsz = 0.0245436926 //0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 40001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 0.0
diff --git a/config/templates/sinktest/256res_rotating_highvisc.conf b/config/templates/sinktest/256res_rotating_highvisc.conf
new file mode 100644
index 0000000..484bb82
--- /dev/null
+++ b/config/templates/sinktest/256res_rotating_highvisc.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 256
+AC_ny = 256 
+AC_nz = 256
+
+AC_dsx = 0.0245436926 //0.04908738521
+AC_dsy = 0.0245436926 //0.04908738521
+AC_dsz = 0.0245436926 //0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 40001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 2e10
diff --git a/config/templates/sinktest/256res_rotating_lowvisc.conf b/config/templates/sinktest/256res_rotating_lowvisc.conf
new file mode 100644
index 0000000..2dc7082
--- /dev/null
+++ b/config/templates/sinktest/256res_rotating_lowvisc.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 256
+AC_ny = 256 
+AC_nz = 256
+
+AC_dsx = 0.0245436926 //0.04908738521
+AC_dsy = 0.0245436926 //0.04908738521
+AC_dsz = 0.0245436926 //0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 40001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 2.5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 2e10
diff --git a/config/templates/sinktest/512res_non-rotating_highvisc.conf b/config/templates/sinktest/512res_non-rotating_highvisc.conf
new file mode 100644
index 0000000..f58c76b
--- /dev/null
+++ b/config/templates/sinktest/512res_non-rotating_highvisc.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 512
+AC_ny = 512 
+AC_nz = 512
+
+AC_dsx = 0.0122718463  //0.0245436926 //0.04908738521
+AC_dsy = 0.0122718463  //0.0245436926 //0.04908738521
+AC_dsz = 0.0122718463  //0.0245436926 //0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 80001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 0.0
diff --git a/config/templates/sinktest/512res_rotating_highvisc.conf b/config/templates/sinktest/512res_rotating_highvisc.conf
new file mode 100644
index 0000000..c3a9420
--- /dev/null
+++ b/config/templates/sinktest/512res_rotating_highvisc.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 512
+AC_ny = 512 
+AC_nz = 512
+
+AC_dsx = 0.0122718463  //0.0245436926 //0.04908738521
+AC_dsy = 0.0122718463  //0.0245436926 //0.04908738521
+AC_dsz = 0.0122718463  //0.0245436926 //0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 80001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 2e10
diff --git a/config/templates/sinktest/512res_rotating_lowvisc.conf b/config/templates/sinktest/512res_rotating_lowvisc.conf
new file mode 100644
index 0000000..7b60afe
--- /dev/null
+++ b/config/templates/sinktest/512res_rotating_lowvisc.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 512
+AC_ny = 512 
+AC_nz = 512
+
+AC_dsx = 0.0122718463  //0.0245436926 //0.04908738521
+AC_dsy = 0.0122718463  //0.0245436926 //0.04908738521
+AC_dsz = 0.0122718463  //0.0245436926 //0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 80001
+AC_save_steps = 50
+AC_bin_steps = 1000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 1.25e-4
+AC_cs_sound = 1.0
+AC_zeta = 0.0
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-4 //1e-5 
+AC_soft = 0.36
+AC_accretion_range = 0.25
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 2e10
diff --git a/config/templates/six_hour_stable.conf b/config/templates/six_hour_stable.conf
new file mode 100644
index 0000000..2a7791b
--- /dev/null
+++ b/config/templates/six_hour_stable.conf
@@ -0,0 +1,75 @@
+
+
+/*
+ * =============================================================================
+ * "Compile-time" params
+ * =============================================================================
+ */
+AC_nx = 128
+AC_ny = 128
+AC_nz = 128
+
+AC_dsx = 0.04908738521
+AC_dsy = 0.04908738521
+AC_dsz = 0.04908738521
+
+/*
+ * =============================================================================
+ * Run-time params
+ * =============================================================================
+ */
+AC_max_steps = 1800001
+AC_save_steps = 100
+AC_bin_steps = 2000
+AC_bin_save_t = 1e666
+
+// Hydro
+AC_cdt = 0.4
+AC_cdtv = 0.3
+AC_cdts = 1.0
+AC_nu_visc  = 5e-3
+AC_cs_sound = 1.0
+AC_zeta = 0.01
+
+// Magnetic
+AC_eta = 5e-3
+AC_mu0 = 1.4
+AC_chi = 0.0001
+
+// Forcing
+AC_relhel = 0.0
+AC_forcing_magnitude = 1e-20
+AC_kmin              = 0.8
+AC_kmax              = 1.2
+
+
+// Entropy
+AC_cp_sound = 1.0
+AC_gamma = 0.5
+AC_lnT0 = 1.2
+AC_lnrho0 = 1.3
+
+// Sink Particle
+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-10
+AC_soft = 0.36
+AC_accretion_range = 0.15
+
+// Physical properties of the domain
+// Typical 1km/s velocity
+AC_unit_velocity = 1e5
+// using density estimate of 100 H2 molecules per cm^3.
+AC_unit_density = 3e-20
+// using 100,000 A.U.
+AC_unit_length = 1.5e17
+/*
+ * =============================================================================
+ * Initial conditions
+ * =============================================================================
+ */
+AC_ampl_lnrho = 0.0
+AC_ampl_uu = 0.0
diff --git a/doc/manual/sink_particle.md b/doc/manual/sink_particle.md
index eb86610..0b7f11b 100644
--- a/doc/manual/sink_particle.md
+++ b/doc/manual/sink_particle.md
@@ -81,3 +81,103 @@ Make is possible for the particle to accrete momentum in addition to mass, and t
 Create sink particles dynamically and allow the presence of multiple sinks. 
 
 
+### Suggestion writen by JP: 
+
+Add to ```acc/mhd_solver/stencil_defines.h```:
+
+```
+    // Scalar mass
+    #define AC_FOR_USER_REAL_PARAM_TYPES(FUNC) \
+                ...
+                ...
+                FUNC(AC_sink_particle_mass),
+
+    // Vector position
+    // NOTE: This makes an AcReal3 constant parameter
+    // This is a completely new type that has not yet been
+    // tested. Accessible from the DSL with
+    // DCONST_REAL3(AC_sink_particle_pos)
+    #define AC_FOR_USER_REAL3_PARAM_TYPES(FUNC) \
+                ...
+                ...
+                FUNC(AC_sink_particle_pos),
+```
+Currently we do not do it like that as AcReal3 constant parameter. However, would be good thing to implement eventually. 
+```
+
+    // Vertex buffer for accretion
+    #define AC_FOR_VTXBUF_HANDLES(FUNC) \
+                ...
+                FUNC(VTXBUF_ACCRETION),
+```
+
+Add to ```acc/mhd_solver/stencil_process.sps```:
+
+```
+
+    Scalar
+    your_accretion_function(int3 globalVertexIdx)
+    {
+        Scalar mass = DCONST_REAL(AC_sink_particle_mass);
+        Vector pos0 = DCONST_REAL3(AC_sink_particle_pos);
+        Vector pos1 = (Vector){ (globalVertexIdx.x - nx_min) * dsx), ...};
+        
+        return *accretion from the current cell at globalVertexIdx*
+    }
+    
+```
+This step will require that we calculate the acctetion as **mass** in a grid cell volume. This is to eventually take an advantage of nonunifrorm grid. OIOn that way we do not need to worry about it in the collective operation stage. 
+```
+
+    // This should have a global scope
+    out Scalar out_accretion = VTXBUF_ACCRETION;
+
+    Kernel void
+    solve(Scalar dt) {
+
+        ...
+        ...
+
+        out_accretion = your_accretion_function(globalVertexIdx);            
+    }
+```
+Important to note that we need to take into account the reduction of density with other equations of motion. Or apply analogous style to forcing.  
+
+Create new file ```src/standalone/model/host_accretion.cc```:
+
+```
+    #include "astaroth.h"
+
+    void
+    updateAccretion(AcMeshInfo* info)
+    {
+        AcReal previous_mass = info->real_params[AC_sink_particle_mass];
+        AcReal accreted_mass = acReduceScal(RTYPE_SUM, VTXBUF_ACCRETION);
+        // Note: RTYPE_SUM does not yet exist (but is easy to add)
+       
+        AcReal mass = previous_mass + accreted_mass;
+        info->real_params[AC_sink_particle_mass] = mass; // Save for the next iteration
+        acLoadDeviceConstant(AC_sink_particle_mass, mass); // Load to the GPUs
+    }
+```
+We assume that acReduceScal will sum up the mass in the accretion buffer. 
+
+Call from```src/standalone/simulation.cc```:
+
+```    
+    #include "model/host_accretion.h"
+    
+    int
+    run_simulation(void)
+    {
+        ...
+
+        while (simulation_running) {
+            ...
+            ...
+
+            updateAccretion(&mesh_info);
+        }
+    }
+```
+
diff --git a/src/standalone/config_loader.cc b/src/standalone/config_loader.cc
index 36fb83a..98098ff 100644
--- a/src/standalone/config_loader.cc
+++ b/src/standalone/config_loader.cc
@@ -125,22 +125,25 @@ update_config(AcMeshInfo* config)
                                        config->real_params[AC_gamma];
 
     AcReal G_CONST_CGS = AcReal(
-        6.674e-8);                   // g/cm3/s GGS definition //TODO define in a separate module
+        6.674e-8); // cm^3/(g*s^2) GGS definition //TODO define in a separate module
     AcReal M_sun = AcReal(1.989e33); // g solar mass
 
-    config->real_params[AC_M_star] = config->real_params[AC_M_star] * M_sun /
-                                     ((config->real_params[AC_unit_length] *
-                                       config->real_params[AC_unit_length] *
-                                       config->real_params[AC_unit_length]) *
-                                      config->real_params[AC_unit_density]);
+    config->real_params[AC_unit_mass] = (config->real_params[AC_unit_length] *
+                                         config->real_params[AC_unit_length] *
+                                         config->real_params[AC_unit_length]) *
+                                        config->real_params[AC_unit_density];
 
-    config->real_params[AC_G_CONST] = G_CONST_CGS / ((config->real_params[AC_unit_velocity] *
+    config->real_params[AC_M_sink] = config->real_params[AC_M_sink_Msun] * M_sun /
+                                     config->real_params[AC_unit_mass];
+    config->real_params[AC_M_sink_init] = config->real_params[AC_M_sink_Msun] * M_sun /
+                                          config->real_params[AC_unit_mass];
+
+    config->real_params[AC_G_const] = G_CONST_CGS / ((config->real_params[AC_unit_velocity] *
                                                       config->real_params[AC_unit_velocity]) /
                                                      (config->real_params[AC_unit_density] *
+                                                      config->real_params[AC_unit_length] *
                                                       config->real_params[AC_unit_length]));
 
-    config->real_params[AC_GM_star] = config->real_params[AC_M_star] *
-                                      config->real_params[AC_G_CONST];
     config->real_params[AC_sq2GM_star] = AcReal(sqrt(AcReal(2) * config->real_params[AC_GM_star]));
 
 #if VERBOSE_PRINTING // Defined in astaroth.h
diff --git a/src/standalone/model/host_memory.cc b/src/standalone/model/host_memory.cc
index 5a68d9b..f6a54cb 100644
--- a/src/standalone/model/host_memory.cc
+++ b/src/standalone/model/host_memory.cc
@@ -71,7 +71,7 @@ acmesh_create(const AcMeshInfo& mesh_info)
     return mesh;
 }
 
-static void
+void
 vertex_buffer_set(const VertexBufferHandle& key, const AcReal& val, AcMesh* mesh)
 {
     const int n = acVertexBufferSize(mesh->info);
@@ -207,6 +207,83 @@ inflow_vedge(AcMesh* mesh)
     }
 }
 
+// This is the initial condition type for the infalling vedge in the pseudodisk
+// model.
+void
+simple_uniform_core(AcMesh* mesh)
+{
+    const int mx = mesh->info.int_params[AC_mx];
+    const int my = mesh->info.int_params[AC_my];
+    const int mz = mesh->info.int_params[AC_mz];
+
+    const double DX = mesh->info.real_params[AC_dsx];
+    const double DY = mesh->info.real_params[AC_dsy];
+    const double DZ = mesh->info.real_params[AC_dsz];
+
+    const double ampl_lnrho = mesh->info.real_params[AC_ampl_lnrho];
+
+    const double xorig   = mesh->info.real_params[AC_xorig];
+    const double yorig   = mesh->info.real_params[AC_yorig];
+    const double zorig   = mesh->info.real_params[AC_zorig];
+
+    const double G_const = mesh->info.real_params[AC_G_const]; 
+    const double M_sink_init = mesh->info.real_params[AC_M_sink_init];
+    const double cs2_sound = mesh->info.real_params[AC_cs2_sound];
+
+    const double RR_inner_bound = mesh->info.real_params[AC_soft]/AcReal(2.0);
+    const double core_coeff   = (exp(ampl_lnrho)  * cs2_sound) / (double(4.0)*M_PI * G_const);
+
+    double xx, yy, zz, RR;
+    double core_profile;
+
+    //TEMPORARY TEST INPUT PARAMETERS
+    const double core_radius = DX*32.0;
+    const double trans = DX*12.0;
+    //const double epsilon = DX*2.0;
+    const double vel_scale = mesh->info.real_params[AC_ampl_uu];
+    double abso_vel;
+
+    RR = 1.0;
+    printf("%e %e %e \n", RR, trans, core_radius);
+     
+
+    for (int k = 0; k < mz; k++) {
+        for (int j = 0; j < my; j++) {
+            for (int i = 0; i < mx; i++) {
+                int idx = i + j*mx + k*mx*my;
+                xx     = DX * double(i) - xorig;
+                yy     = DY * double(j) - yorig;
+                zz     = DZ * double(k) - zorig;
+
+                RR     = sqrt(xx*xx + yy*yy + zz*zz);
+                     
+                if (RR >= RR_inner_bound) {
+                    abso_vel = vel_scale * sqrt(2.0 * G_const 
+                                                    * M_sink_init / RR);
+                    core_profile = pow(RR, -2.0);   //double(1.0);
+                } else {
+                    abso_vel = vel_scale * sqrt(2.0 * AC_G_const 
+                                                    * AC_M_sink_init / RR_inner_bound);
+                    core_profile = pow(RR_inner_bound, -2.0);   //double(1.0);
+                }
+
+               if (RR <= sqrt(DX*DX + DY*DY + DZ*DZ)) {
+                    abso_vel = 0.0;
+                    RR = 1.0;
+               }
+
+
+                mesh->vertex_buffer[VTXBUF_LNRHO][idx] = AcReal(log(core_coeff*core_profile));
+                mesh->vertex_buffer[VTXBUF_UUX][idx]   = AcReal(-abso_vel * (yy / RR));
+                mesh->vertex_buffer[VTXBUF_UUY][idx]   = AcReal( abso_vel * (xx / RR));
+                mesh->vertex_buffer[VTXBUF_UUZ][idx]   = AcReal(0.0);
+
+            }
+        }
+    }
+}
+
+
 // This is the initial condition type for the infalling vedge in the pseudodisk
 // model.
 void
@@ -564,6 +641,7 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
     }
     case INIT_TYPE_GAUSSIAN_RADIAL_EXPL:
         acmesh_clear(mesh);
+        vertex_buffer_set(VTXBUF_LNRHO, mesh->info.real_params[AC_ampl_lnrho], mesh);
         // acmesh_init_to(INIT_TYPE_RANDOM, mesh);
         gaussian_radial_explosion(mesh);
 
@@ -581,6 +659,10 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
             }
         }
         break;
+    case INIT_TYPE_SIMPLE_CORE:
+        acmesh_clear(mesh);
+        simple_uniform_core(mesh);
+        break;
     case INIT_TYPE_VEDGE:
         acmesh_clear(mesh);
         inflow_vedge_freefall(mesh);
diff --git a/src/standalone/model/host_memory.h b/src/standalone/model/host_memory.h
index e2958cd..b4f1a77 100644
--- a/src/standalone/model/host_memory.h
+++ b/src/standalone/model/host_memory.h
@@ -34,6 +34,7 @@
         FUNC(INIT_TYPE_XWAVE), \
         FUNC(INIT_TYPE_GAUSSIAN_RADIAL_EXPL), \
         FUNC(INIT_TYPE_ABC_FLOW) , \
+        FUNC(INIT_TYPE_SIMPLE_CORE), \
         FUNC(INIT_TYPE_VEDGE), \
         FUNC(INIT_TYPE_VEDGEX), \
         FUNC(INIT_TYPE_RAYLEIGH_TAYLOR), \
@@ -48,6 +49,8 @@ extern const char* init_type_names[]; // Defined in host_memory.cc
 
 AcMesh* acmesh_create(const AcMeshInfo& mesh_info);
 
+void vertex_buffer_set(const VertexBufferHandle& key, const AcReal& val, AcMesh* mesh);
+
 void acmesh_clear(AcMesh* mesh);
 
 void acmesh_init_to(const InitType& type, AcMesh* mesh);
diff --git a/src/standalone/renderer.cc b/src/standalone/renderer.cc
index 2e8945f..be36a37 100644
--- a/src/standalone/renderer.cc
+++ b/src/standalone/renderer.cc
@@ -41,6 +41,10 @@
 #include "src/core/math_utils.h"
 #include "timer_hires.h"
 
+// NEED TO BE DEFINED HERE. IS NOT NOTICED BY compile_acc call.
+#define LFORCING (1)
+#define LSINK (0)
+
 // Window
 SDL_Renderer* renderer      = NULL;
 static SDL_Window* window   = NULL;
@@ -293,6 +297,7 @@ renderer_quit(void)
 }
 
 static int init_type = INIT_TYPE_GAUSSIAN_RADIAL_EXPL;
+// static int init_type = INIT_TYPE_SIMPLE_CORE;
 
 static bool
 running(AcMesh* mesh)
@@ -359,6 +364,10 @@ run_renderer(void)
     AcMesh* mesh = acmesh_create(mesh_info);
     acmesh_init_to(InitType(init_type), mesh);
 
+#if LSINK
+    vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh);
+#endif
+
     acInit(mesh_info);
     acLoad(*mesh);
 
@@ -375,6 +384,9 @@ run_renderer(void)
     int steps                      = 0;
     k_slice                        = mesh->info.int_params[AC_mz] / 2;
     k_slice_max                    = mesh->info.int_params[AC_mz];
+#if LSINK
+    AcReal accreted_mass = 0.0;
+#endif
     while (running(mesh)) {
 
         /* Input */
@@ -382,15 +394,25 @@ run_renderer(void)
         timer_reset(&io_timer);
 
 /* Step the simulation */
-#if 1
-        const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
-        const AcReal dt   = host_timestep(umax, mesh_info);
+#if LSINK
+        const AcReal sum_mass = acReduceScal(RTYPE_SUM, VTXBUF_ACCRETION);
+        accreted_mass         = accreted_mass + sum_mass;
+        AcReal sink_mass      = mesh_info.real_params[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);
+        vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh);
+#endif
 
 #if LFORCING
         const ForcingParams forcing_params = generateForcingParams(mesh->info);
         loadForcingParamsToDevice(forcing_params);
 #endif
 
+#if 1
+        const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
+        const AcReal dt   = host_timestep(umax, mesh_info);
+
         acIntegrate(dt);
 #else
         ModelMesh* model_mesh = modelmesh_create(mesh->info);
diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc
index 51a9a76..bdf0b18 100644
--- a/src/standalone/simulation.cc
+++ b/src/standalone/simulation.cc
@@ -40,6 +40,10 @@
 #include <sys/stat.h>
 #include <sys/types.h>
 
+// NEED TO BE DEFINED HERE. IS NOT NOTICED BY compile_acc call.
+#define LFORCING (1)
+#define LSINK (0)
+
 // Write all setting info into a separate ascii file. This is done to guarantee
 // that we have the data specifi information in the thing, even though in
 // principle these things are in the astaroth.conf.
@@ -56,37 +60,76 @@ write_mesh_info(const AcMeshInfo* config)
 
     infotxt = fopen("mesh_info.list", "w");
 
+    /*
+    // JP: this could be done shorter and with smaller chance for errors with the following
+    // (modified from acPrintMeshInfo() in astaroth.cu)
+
+    for (int i = 0; i < NUM_INT_PARAMS; ++i)
+        fprintf(infotxt, "int %s: %d\n", intparam_names[i], config.int_params[i]);
+
+    for (int i = 0; i < NUM_INT3_PARAMS; ++i)
+        fprintf(infotxt, "int3 %s: (%d, %d, %d)\n", int3param_names[i], config.int3_params[i].x,
+                                                                        config.int3_params[i].y,
+                                                                        config.int3_params[i].z);
+
+    for (int i = 0; i < NUM_REAL_PARAMS; ++i)
+        fprintf(infotxt, "real %s: %g\n", realparam_names[i], double(config.real_params[i]));
+
+    for (int i = 0; i < NUM_REAL3_PARAMS; ++i)
+        fprintf(infotxt, "real %s: (%g, %g, %g)\n", real3param_names[i],
+                                                    double(config.real3_params[i].x),
+                                                    double(config.real3_params[i].y),
+                                                    double(config.real3_params[i].z));
+    */
+
     // Total grid dimensions
-    fprintf(infotxt, "int  AC_mx        %i \n", config->int_params[AC_mx]);
-    fprintf(infotxt, "int  AC_my        %i \n", config->int_params[AC_my]);
-    fprintf(infotxt, "int  AC_mz        %i \n", config->int_params[AC_mz]);
+    fprintf(infotxt, "int  AC_mx         %i \n", config->int_params[AC_mx]);
+    fprintf(infotxt, "int  AC_my         %i \n", config->int_params[AC_my]);
+    fprintf(infotxt, "int  AC_mz         %i \n", config->int_params[AC_mz]);
 
     // Bounds for the computational domain, i.e. nx_min <= i < nx_max
-    fprintf(infotxt, "int  AC_nx_min    %i \n", config->int_params[AC_nx_min]);
-    fprintf(infotxt, "int  AC_nx_max    %i \n", config->int_params[AC_nx_max]);
-    fprintf(infotxt, "int  AC_ny_min    %i \n", config->int_params[AC_ny_min]);
-    fprintf(infotxt, "int  AC_ny_max    %i \n", config->int_params[AC_ny_max]);
-    fprintf(infotxt, "int  AC_nz_min    %i \n", config->int_params[AC_nz_min]);
-    fprintf(infotxt, "int  AC_nz_max    %i \n", config->int_params[AC_nz_max]);
+    fprintf(infotxt, "int  AC_nx_min     %i \n", config->int_params[AC_nx_min]);
+    fprintf(infotxt, "int  AC_nx_max     %i \n", config->int_params[AC_nx_max]);
+    fprintf(infotxt, "int  AC_ny_min     %i \n", config->int_params[AC_ny_min]);
+    fprintf(infotxt, "int  AC_ny_max     %i \n", config->int_params[AC_ny_max]);
+    fprintf(infotxt, "int  AC_nz_min     %i \n", config->int_params[AC_nz_min]);
+    fprintf(infotxt, "int  AC_nz_max     %i \n", config->int_params[AC_nz_max]);
 
     // Spacing
-    fprintf(infotxt, "real AC_dsx       %e \n", (double)config->real_params[AC_dsx]);
-    fprintf(infotxt, "real AC_dsy       %e \n", (double)config->real_params[AC_dsy]);
-    fprintf(infotxt, "real AC_dsz       %e \n", (double)config->real_params[AC_dsz]);
-    fprintf(infotxt, "real AC_inv_dsx   %e \n", (double)config->real_params[AC_inv_dsx]);
-    fprintf(infotxt, "real AC_inv_dsy   %e \n", (double)config->real_params[AC_inv_dsy]);
-    fprintf(infotxt, "real AC_inv_dsz   %e \n", (double)config->real_params[AC_inv_dsz]);
-    fprintf(infotxt, "real AC_dsmin     %e \n", (double)config->real_params[AC_dsmin]);
+    fprintf(infotxt, "real AC_dsx        %e \n", (double)config->real_params[AC_dsx]);
+    fprintf(infotxt, "real AC_dsy        %e \n", (double)config->real_params[AC_dsy]);
+    fprintf(infotxt, "real AC_dsz        %e \n", (double)config->real_params[AC_dsz]);
+    fprintf(infotxt, "real AC_inv_dsx    %e \n", (double)config->real_params[AC_inv_dsx]);
+    fprintf(infotxt, "real AC_inv_dsy    %e \n", (double)config->real_params[AC_inv_dsy]);
+    fprintf(infotxt, "real AC_inv_dsz    %e \n", (double)config->real_params[AC_inv_dsz]);
+    fprintf(infotxt, "real AC_dsmin      %e \n", (double)config->real_params[AC_dsmin]);
 
     /* Additional helper params */
     // Int helpers
-    fprintf(infotxt, "int  AC_mxy       %i \n", config->int_params[AC_mxy]);
-    fprintf(infotxt, "int  AC_nxy       %i \n", config->int_params[AC_nxy]);
-    fprintf(infotxt, "int  AC_nxyz      %i \n", config->int_params[AC_nxyz]);
+    fprintf(infotxt, "int  AC_mxy        %i \n", config->int_params[AC_mxy]);
+    fprintf(infotxt, "int  AC_nxy        %i \n", config->int_params[AC_nxy]);
+    fprintf(infotxt, "int  AC_nxyz       %i \n", config->int_params[AC_nxyz]);
 
     // Real helpers
-    fprintf(infotxt, "real AC_cs2_sound %e \n", (double)config->real_params[AC_cs2_sound]);
-    fprintf(infotxt, "real AC_cv_sound  %e \n", (double)config->real_params[AC_cv_sound]);
+    fprintf(infotxt, "real AC_cs2_sound  %e \n", (double)config->real_params[AC_cs2_sound]);
+    fprintf(infotxt, "real AC_cv_sound   %e \n", (double)config->real_params[AC_cv_sound]);
+
+    // Physical units
+    fprintf(infotxt, "real AC_unit_density  %e \n", (double)config->real_params[AC_unit_density]);
+    fprintf(infotxt, "real AC_unit_velocity %e \n", (double)config->real_params[AC_unit_velocity]);
+    fprintf(infotxt, "real AC_unit_mass     %e \n", (double)config->real_params[AC_unit_mass]);
+    fprintf(infotxt, "real AC_unit_length   %e \n", (double)config->real_params[AC_unit_length]);
+
+    // Here I'm still trying to copy the structure of the code above, and see if this will work for
+    // sink particle. I haven't fully undertand what these lines do but I'll read up on them soon.
+    // This is still yet experimental.
+    // Sink particle
+    fprintf(infotxt, "real AC_sink_pos_x %e \n", (double)config->real_params[AC_sink_pos_x]);
+    fprintf(infotxt, "real AC_sink_pos_y %e \n", (double)config->real_params[AC_sink_pos_y]);
+    fprintf(infotxt, "real AC_sink_pos_z %e \n", (double)config->real_params[AC_sink_pos_z]);
+    fprintf(infotxt, "real AC_M_sink     %e \n", (double)config->real_params[AC_M_sink]);
+    fprintf(infotxt, "real AC_soft     %e \n", (double)config->real_params[AC_soft]);
+    fprintf(infotxt, "real AC_G_const     %e \n", (double)config->real_params[AC_G_const]);
 
     fclose(infotxt);
 }
@@ -135,7 +178,8 @@ save_mesh(const AcMesh& save_mesh, const int step, const AcReal t_step)
 // This function prints out the diagnostic values to std.out and also saves and
 // appends an ascii file to contain all the result.
 static inline void
-print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* diag_file)
+print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* diag_file,
+                  const AcReal sink_mass, const AcReal accreted_mass)
 {
 
     AcReal buf_rms, buf_max, buf_min;
@@ -165,6 +209,10 @@ print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* di
         fprintf(diag_file, "%e %e %e ", double(buf_min), double(buf_rms), double(buf_max));
     }
 
+    if ((sink_mass >= AcReal(0.0)) || (accreted_mass >= AcReal(0.0))) {
+        fprintf(diag_file, "%e %e ", double(sink_mass), double(accreted_mass));
+    }
+
     fprintf(diag_file, "\n");
 }
 
@@ -184,6 +232,11 @@ run_simulation(void)
     AcMesh* mesh = acmesh_create(mesh_info);
     // TODO: This need to be possible to define in astaroth.conf
     acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, mesh);
+    // acmesh_init_to(INIT_TYPE_SIMPLE_CORE, mesh); //Initial condition for a collapse test
+
+#if LSINK
+    vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh);
+#endif
 
     acInit(mesh_info);
     acLoad(*mesh);
@@ -199,11 +252,17 @@ run_simulation(void)
         fprintf(diag_file, "%s_min  %s_rms  %s_max  ", vtxbuf_names[i], vtxbuf_names[i],
                 vtxbuf_names[i]);
     }
-
+#if LSINK
+    fprintf(diag_file, "sink_mass  accreted_mass  ");
+#endif
     fprintf(diag_file, "\n");
 
     write_mesh_info(&mesh_info);
-    print_diagnostics(0, AcReal(.0), t_step, diag_file);
+#if LSINK
+    print_diagnostics(0, AcReal(.0), t_step, diag_file, mesh_info.real_params[AC_M_sink_init], 0.0);
+#else
+    print_diagnostics(0, AcReal(.0), t_step, diag_file, -1.0, -1.0);
+#endif
 
     acBoundcondStep();
     acStore(mesh);
@@ -228,35 +287,52 @@ run_simulation(void)
     //  acUpdate_sink_particle() will do the similar trick to the device.
 
     /* Step the simulation */
+    AcReal accreted_mass = 0.0;
+    AcReal sink_mass     = 0.0;
     for (int i = 1; i < max_steps; ++i) {
         const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
         const AcReal dt   = host_timestep(umax, mesh_info);
 
+#if LSINK
+
+        const AcReal sum_mass = acReduceScal(RTYPE_SUM, VTXBUF_ACCRETION);
+        accreted_mass         = accreted_mass + sum_mass;
+        sink_mass             = 0.0;
+        sink_mass             = mesh_info.real_params[AC_M_sink_init] + accreted_mass;
+        acLoadDeviceConstant(AC_M_sink, sink_mass);
+        vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh);
+
+        int on_off_switch;
+        if (i < 1) {
+            on_off_switch = 0; // accretion is off till certain amount of steps.
+        }
+        else {
+            on_off_switch = 1;
+        }
+        acLoadDeviceConstant(AC_switch_accretion, on_off_switch);
+
+        // MV: Old TODOs to remind of eventual future directions.
+        // TODO_SINK acUpdate_sink_particle()
+        //  3. Velocity of the particle)
+        // TODO_SINK acAdvect_sink_particle()
+        //  1. Calculate the equation of motion for the sink particle.
+        //  NOTE: Might require embedding with acIntegrate(dt).
+        // TODO_SINK acAccrete_sink_particle()
+        //  2. Transfer momentum into sink particle
+        //  (OPTIONAL: Affection the motion of the particle)
+        //  NOTE: Might require embedding with acIntegrate(dt).
+        //  This is the hardest part. Please see Lee et al. ApJ 783 (2014) for reference.
+#else
+        accreted_mass = -1.0;
+        sink_mass     = -1.0;
+#endif
+
 #if LFORCING
         const ForcingParams forcing_params = generateForcingParams(mesh_info);
         loadForcingParamsToDevice(forcing_params);
 #endif
 
-        // TODO_SINK acUpdate_sink_particle()
-        //  Update properties of the sing particle for acIntegrate(). Essentially:
-        //  1. Location of the particle
-        //  2. Mass of the particle
-        //  (3. Velocity of the particle)
-        //  These can be used for calculating he gravitational field.
-
         acIntegrate(dt);
-        // TODO_SINK acAdvect_sink_particle()
-        //  THIS IS OPTIONAL. We will start from unmoving particle.
-        //  1. Calculate the equation of motion for the sink particle.
-        //  NOTE: Might require embedding with acIntegrate(dt).
-
-        // TODO_SINK acAccrete_sink_particle()
-        //  Calculate accretion of the sink particle from the surrounding medium
-        //  1. Transfer density into sink particle mass
-        //  2. Transfer momentum into sink particle
-        //  (OPTIONAL: Affection the motion of the particle)
-        //  NOTE: Might require embedding with acIntegrate(dt).
-        //  This is the hardest part. Please see Lee et al. ApJ 783 (2014) for reference.
 
         t_step += dt;
 
@@ -269,8 +345,11 @@ run_simulation(void)
                 timeseries.ts.
             */
 
-            print_diagnostics(i, dt, t_step, diag_file);
-
+            print_diagnostics(i, dt, t_step, diag_file, sink_mass, accreted_mass);
+#if LSINK
+            printf("sink mass is: %.15e \n", double(sink_mass));
+            printf("accreted mass is: %.15e \n", double(accreted_mass));
+#endif
             /*
                 We would also might want an XY-average calculating funtion,
                 which can be very useful when observing behaviour of turbulent