From 632819764eb964650d2da02bf6d582829c4923f3 Mon Sep 17 00:00:00 2001 From: JackHsu Date: Tue, 6 Aug 2019 11:47:55 +0800 Subject: [PATCH] Some intuitions about accretion process. Also defined as a module. --- acc/mhd_solver/stencil_defines.h | 1 + acc/mhd_solver/stencil_process.sps | 23 +++++++++++++++++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/acc/mhd_solver/stencil_defines.h b/acc/mhd_solver/stencil_defines.h index bd6e43c..d37eeba 100644 --- a/acc/mhd_solver/stencil_defines.h +++ b/acc/mhd_solver/stencil_defines.h @@ -33,6 +33,7 @@ #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 6ac6e4f..ea82ca5 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -56,7 +56,6 @@ continuity(in Vector uu, in Scalar lnrho) { } #if LSINK -//first attempt to do a self-containing LSINK module Vector sink_gravity(int3 globalVertexIdx){ Vector force_gravity; @@ -65,8 +64,8 @@ sink_gravity(int3 globalVertexIdx){ (globalVertexIdx.z - nz_min) * dsz}; const Scalar sink_mass = DCONST_REAL(AC_M_sink); const Vector sink_pos = (Vector){DCONST_REAL(AC_sink_pos_x), - DCONST_REAL(AC_sink_pos_y), - DCONST_REAL(AC_sink_pos_z)}; + DCONST_REAL(AC_sink_pos_y), + DCONST_REAL(AC_sink_pos_z)}; const Scalar distance = length(grid_pos - sink_pos); const Scalar soft = DCONST_REAL(AC_soft); const Scalar gravity_magnitude = (AC_G_const * sink_mass) / pow(((distance * distance) + soft*soft), 1.5); @@ -77,7 +76,23 @@ sink_gravity(int3 globalVertexIdx){ return force_gravity; } #endif - +// Now this is more of a suedo-code, just write out some intuitions. +#if LACCRETION +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; +} #if LENTROPY Vector momentum(int3 globalVertexIdx, in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {