Merged in sink_20190723 (pull request #8)

Sink 20190723
2019-09-17 15:14:34 +00:00
parent 55e4357d77 d2057b8f5a
commit 4b6ce26f8a
25 changed files with 1564 additions and 121 deletions
--- 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
+
--- 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
 }
--- 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")
--- 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()
--- 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,
--- 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)


--- 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
--- a/config/templates/core_collapse.conf
+++ 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
--- a/config/templates/sink_aftersummer/128_basic_rotating.conf
+++ 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
--- a/config/templates/sinktest.conf
+++ 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
--- a/config/templates/sinktest/128res_non-rotating_highvisc.conf
+++ 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
--- a/config/templates/sinktest/128res_rotating_highvisc.conf
+++ 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
--- a/config/templates/sinktest/256res_non-rotating_highvisc.conf
+++ 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
--- a/config/templates/sinktest/256res_rotating_highvisc.conf
+++ 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
--- a/config/templates/sinktest/256res_rotating_lowvisc.conf
+++ 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
--- a/config/templates/sinktest/512res_non-rotating_highvisc.conf
+++ 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
--- a/config/templates/sinktest/512res_rotating_highvisc.conf
+++ 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
--- a/config/templates/sinktest/512res_rotating_lowvisc.conf
+++ 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
--- a/config/templates/six_hour_stable.conf
+++ 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
--- 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);
+        }
+    }
+```
+
--- 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
--- 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);
--- 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);
--- 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);
--- 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