cwpearson/astaroth
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added update_accretion_buffer.

Browse Source
This commit is contained in:
JackHsu
2019-08-08 15:53:41 +08:00
parent c1b99b2b37
commit 89128af44b
3 changed files with 24 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
acc/mhd_solver/stencil_defines.h
View File

@@ -80,6 +80,7 @@
FUNC(AC_M_sink),\ FUNC(AC_M_sink),\
FUNC(AC_M_sink_Msun),\ FUNC(AC_M_sink_Msun),\
FUNC(AC_soft),\ FUNC(AC_soft),\
FUNC(AC_accretion_range),\
/* Run params */\ /* Run params */\
FUNC(AC_cdt), \ FUNC(AC_cdt), \
FUNC(AC_cdtv), \ FUNC(AC_cdtv), \

43
acc/mhd_solver/stencil_process.sps
View File

@@ -82,14 +82,14 @@ sink_gravity(int3 globalVertexIdx){
#if LSINK #if LSINK
Scalar Scalar
truelove_density(in Scalar lnrho){ truelove_density(in Scalar lnrho){
const Scalar rho = exp(lnrho); const Scalar rho = exp(value(lnrho));
const Scalar Jeans_length_squared = (M_PI * cs2_sound) / (AC_G_const * rho); 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); 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. //TODO: dsx will cancel out, deal with it later for optimization.
Scalar accrection_rho = rho - TJ_rho; Scalar accretion_rho = rho - TJ_rho;
if (accrection_rho < 0){ if (accretion_rho < 0){
accrection_rho = Scalar(0); accretion_rho = Scalar(0);
} }
return accretion_rho; return accretion_rho;
@@ -97,23 +97,24 @@ truelove_density(in Scalar lnrho){
Scalar Scalar
accretion_profile(int3 globalVertexIdx, in Scalar lnrho){ accretion_profile(int3 globalVertexIdx, in Scalar lnrho){
// QUESTION: do I need to define grid_pos, sink_pos and distance again // 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. // 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, const Vector grid_pos = (Vector){(globalVertexIdx.x - nx_min) * dsx,
(globalVertexIdx.y - ny_min) * dsy, (globalVertexIdx.y - ny_min) * dsy,
(globalVertexIdx.z - nz_min) * dsz}; (globalVertexIdx.z - nz_min) * dsz};
const Vector sink_pos = (Vector){DCONST_REAL(AC_sink_pos_x), const Vector sink_pos = (Vector){DCONST_REAL(AC_sink_pos_x),
DCONST_REAL(AC_sink_pos_y), DCONST_REAL(AC_sink_pos_y),
DCONST_REAL(AC_sink_pos_z)}; DCONST_REAL(AC_sink_pos_z)};
const Scalar profile_range = DCONST_REAL(AC_accretion_range); const Scalar profile_range = DCONST_REAL(AC_accretion_range);
const Scalar accretion_distance = length(grid_pos - sink_pos); const Scalar accretion_distance = length(grid_pos - sink_pos);
if ((accretion_distance) <= profile_range){ Scalar accretion_density;
if ((accretion_distance) <= profile_range){
// calculate accretion according to chosen criterion for the grid cell. // calculate accretion according to chosen criterion for the grid cell.
accretion_density = truelove_density(lnrho); accretion_density = truelove_density(lnrho);
} else { } else {
accretion_density = Scalar(0.0); accretion_density = Scalar(0.0);
} }
return accretion_density; return accretion_density;
} }
#endif #endif
@@ -417,7 +418,7 @@ solve(Scalar dt) {
#if LSINK #if LSINK
out_lnrho = log(exp(out_lnrho) - accretion_profile(globalVertexIdx, lnrho)); out_lnrho = log(exp(out_lnrho) - accretion_profile(globalVertexIdx, lnrho));
out_accretion = accretion + (accretion_profile(globalVertexIdx,lnrho) * dsx * dsy * dsz); out_accretion = value(accretion) + (accretion_profile(globalVertexIdx,lnrho) * dsx * dsy * dsz);
#endif #endif
} }

2
config/astaroth.conf
View File

@@ -58,7 +58,7 @@ AC_soft = 0.12
// Accretion Parameters // Accretion Parameters
// profile_range is multiple of dsx // profile_range is multiple of dsx
AC_profile_range = 2.0 AC_accretion_range = 2.0
// Physical properties of the domain // Physical properties of the domain
AC_unit_velocity = 1.0 AC_unit_velocity = 1.0