From aebf4f11aeebc148a8cda12227686be889d3c90a Mon Sep 17 00:00:00 2001
From: JackHsu <ssp201903@asiaa.sinica.edu.tw>
Date: Tue, 6 Aug 2019 16:17:11 +0800
Subject: [PATCH] Added function for truelove-density, and took out the
 unecessary use of switch on-off module for accretion.

---
 acc/mhd_solver/stencil_defines.h   |  1 -
 acc/mhd_solver/stencil_process.sps | 35 ++++++++++++++++++------------
 2 files changed, 21 insertions(+), 15 deletions(-)

diff --git a/acc/mhd_solver/stencil_defines.h b/acc/mhd_solver/stencil_defines.h
index d37eeba..bd6e43c 100644
--- a/acc/mhd_solver/stencil_defines.h
+++ b/acc/mhd_solver/stencil_defines.h
@@ -33,7 +33,6 @@
 #define LFORCING (0)
 #define LUPWD (0)
 #define LSINK (1)
-#define LACCRETION(1)
 
 #define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter
 
diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index ea82ca5..6f469cc 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -76,23 +76,30 @@ sink_gravity(int3 globalVertexIdx){
     return force_gravity;
 }
 #endif
+
+
+
 // Now this is more of a suedo-code, just write out some intuitions. 
-#if LACCRETION
+#if LSINK
 Scalar
-mass_accrection(int3 globalVertexIdx){
-    Scalar mass = (Scalar){()
-    Scalar Jeans_length = pow(((3.14 * cv_sound * cv_sound) / (AC_G_const * lnrho), 0.5));
-    Scalar Truelove_Jeans_lnrho = ((3.14) * Jeans_length * Jeans_length * cv2_sound) / (AC_G_const * dsx * dsy * dsz);    
-    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 sink_mass = DCONST_REAL(AC_M_sink);
-    if (lnrho >= Truelove_Jeans_lnrho)
-       {
-        sink_mass = sink_mass + (Truelove_Jeans_lnrho * (dsx * dsy * dsz))      
-       }
-    return sink_mass;
+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);
+       
+ 
+    Scalar accrection_rho = rho - TJ_rho;
+    if (accrection_rho < 0){
+     accrection_rho = Scalar(0);
+    } 
+       
+    return accretion_rho;
 }
+#endif
+//TODO: basic structure of this part is as follows 
+// update_accretion_buffer() <--> accretion_profile() <--> truelove_density()
+
+
 #if LENTROPY
 Vector
 momentum(int3 globalVertexIdx, in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {