From dac509241fa8f27d267714e679fa9c7fecbabbe6 Mon Sep 17 00:00:00 2001
From: Miikka Vaisala <mvaisala@asiaa.sinica.edu.tw>
Date: Tue, 3 Sep 2019 15:24:58 +0800
Subject: [PATCH] Cleanup, comments and a new template.

---
 acc/mhd_solver/stencil_process.sps            | 15 ++--
 .../sink_aftersummer/128_basic_rotating.conf  | 75 +++++++++++++++++++
 2 files changed, 82 insertions(+), 8 deletions(-)
 create mode 100644 config/templates/sink_aftersummer/128_basic_rotating.conf

diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index ba8bbb7..b888cbc 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -67,6 +67,7 @@ truelove_density(in ScalarField lnrho){
     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,
@@ -82,7 +83,6 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
     Scalar weight;
 
     if (accretion_switch == 1){   
-//    const Scalar weight = exp(-(accretion_distance/profile_range));   
         // Step function weighting 
         if ((accretion_distance) <= profile_range){
             weight = Scalar(1.0);
@@ -90,12 +90,8 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
             weight = Scalar(0.0);
         }
         
-//        const Scalar lnrho_min = Scalar(-10.0);  //TODO Define from astaroth.conf
+        //Truelove criterion is used as a kind of arbitrary density floor. 
         const Scalar lnrho_min = log(truelove_density(lnrho));
-//        const Scalar sink_mass = AC_M_sink;
-//        const Scalar B = Scalar(0.5);
-//        const Scalar k = Scalar(1.5);
-//        const Scalar rate = B * (pow(sink_mass, k) / (AC_dsx * AC_dsy * AC_dsz));
         Scalar rate;     
         if (value(lnrho) > lnrho_min) {
             rate = (exp(value(lnrho)) - exp(lnrho_min)) / dt;
@@ -104,12 +100,12 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
         }
         accretion_density = weight * rate ;
     } else {
-        accretion_density = 0; 
+        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,
@@ -126,6 +122,7 @@ sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) {
     if (accretion_switch == 1){
         Scalar weight;
         // Step function weighting 
+        // MV: This is too aggeressive creating velocity artefacts. 
         if ((accretion_distance) <= profile_range){
         weight = Scalar(1.0);
         } else {
@@ -133,6 +130,8 @@ sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) {
         }
 
         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 { 
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..e6c769d
--- /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 = 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 = 0.0
+AC_ampl_uu = 0.5