From 2f14bb2a30c2fe5043744fc2d443c2a5ce03cc4d Mon Sep 17 00:00:00 2001
From: JackHsu <ssp201903@asiaa.sinica.edu.tw>
Date: Wed, 7 Aug 2019 17:14:26 +0800
Subject: [PATCH] Finished accretion_profile function and started a draft of
 update_accretion_buffer.

---
 acc/mhd_solver/stencil_process.sps | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index 384ecc9..f54cc75 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -85,7 +85,7 @@ truelove_density(in Scalar lnrho){
     const Scalar rho = exp(lnrho);
     const Scalar Jeans_length_squared = (M_PI * cs2_sound) / (AC_G_const * rho);
     const Scalar TJ_rho = ((M_PI) * ((dsx * dsx) / Jeans_length_squared) * cs2_sound) / (AC_G_const * dsx * dsx);
-       
+    //TODO: dsx will cancel out, deal with it later for optimization.      
  
     Scalar accrection_rho = rho - TJ_rho;
     if (accrection_rho < 0){
@@ -96,8 +96,8 @@ truelove_density(in Scalar lnrho){
 }
 
 
-accretion_profile(int3 globalVertexIdx){
-     // ONE QUESTION: do I need to define grid_pos, sink_pos and distance again
+accretion_profile(int3 globalVertexIdx, in Scalar lnrho){
+     // QUESTION: do I need to define grid_pos, sink_pos and distance again
      // if the sink_gravity kernel will also be called once LSINK swtich is on? Seems redundant.
      const Vector grid_pos = (Vector){(globalVertexIdx.x - nx_min) * dsx,
                                       (globalVertexIdx.y - ny_min) * dsy,
@@ -105,17 +105,22 @@ accretion_profile(int3 globalVertexIdx){
      const Vector sink_pos = (Vector){DCONST_REAL(AC_sink_pos_x),
                                       DCONST_REAL(AC_sink_pos_y),
                                       DCONST_REAL(AC_sink_pos_z)};
-     const Scalar profile_range = DCONST_REAL(AC_accretion_range) * dsx;
+     const Scalar profile_range = DCONST_REAL(AC_accretion_range);
      const Scalar accretion_distance = length(grid_pos - sink_pos);   
      if ((accretion_distance) <= profile_range){
         // calculate accretion according to chosen criterion for the grid cell.
-         accretion = accretion_rho * dsx * dsy * dsz;
+         accretion = truelove_density(lnrho);
+     } else {
+         accretion = Scalar(0.0); 
      }
      return accretion;// the returned value is effectively mass, doesn't seem right since MIIKKA said mass will be summed up in vertex buffer. But I'll keep it as it untill we discuss this.
 }
 
 
 update_accretion_buffer(){
+// 1. reduce accretion density from the density field.
+// 2. Add the accretion mass, which is calculated from accretion density times volume of each cell and sum them into accretion buffer.
+
 
 
 
@@ -391,6 +396,11 @@ in Scalar tt = VTXBUF_TEMPERATURE;
 out Scalar out_tt = VTXBUF_TEMPERATURE;
 #endif
 
+#if LSINK
+in Scalar accretion = VTXBUF_ACCRETION;
+out Scalar out_accretion = VTXBUF_ACCRETION;
+#endif
+
 Kernel void
 solve(Scalar dt) {
     out_lnrho = rk3(out_lnrho, lnrho, continuity(uu, lnrho), dt);