From 820132cfe9229de026b2f0c8e7187e82b604b9ec Mon Sep 17 00:00:00 2001 From: JackHsu Date: Mon, 29 Jul 2019 17:53:04 +0800 Subject: [PATCH] I attemped to write a complete LSINK module. Haven't tested how it works but just serve as a save point. --- acc/mhd_solver/stencil_defines.h | 1 + acc/mhd_solver/stencil_process.sps | 25 +++++++++++++++++-------- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/acc/mhd_solver/stencil_defines.h b/acc/mhd_solver/stencil_defines.h index 37cc9e5..c807401 100644 --- a/acc/mhd_solver/stencil_defines.h +++ b/acc/mhd_solver/stencil_defines.h @@ -32,6 +32,7 @@ #define LTEMPERATURE (0) #define LFORCING (0) #define LUPWD (0) +#define LSINK (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 77e55a4..9cf7aa6 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -4,7 +4,7 @@ uniform Scalar nu_visc; uniform Scalar cp_sound; uniform Scalar cv_sound; uniform Scalar mu0; -uniform Scalar eta; +unicalar eta; uniform Scalar gamma; uniform Scalar zeta; @@ -60,13 +60,22 @@ continuity(in Vector uu, in Scalar lnrho) { - divergence(uu); } -//#if SINK -//Here I'm more or less copying the format of the LENTROPY module below -//Vector -//gravity(in Vector xx, in Vector yy, in Vector zz) { - //Here I need to define how gravity works for sink particle, probably based on the distance to the sink particle and its mass. - //Then apply Newton's equation for force, and also remember to consider xyz directions. -//} +#if LSINK + +Vector +gravity() { + vector force_gravity; + const int3 grid_pos = (int3){threadIdx.x + blockIdx.x * blockDim.x + start.x, + threadIdx.y + blockIdx.y * blockDim.y + start.y, + threadIdx.z + blockIdx.z * blockDim.z + start.z}; + const sink_mass = DCONST(AC_M_sink); + const sink_pos = (Vector){DCONST_REAL(AC_sink_pos_x), + DCONST_REAL(AC_sink_pos_y), + DCONST_REAL(AC_sink_pos_z)}; + const distance = dlength_vec(grid_pos - sink_pos); + force_gravity = (6.67e-8 * sink_mass) / (distance * distance); + return force_gravity; +}//first attempt to do a self-containing LSINK module #if LENTROPY Vector