diff --git a/acc/mhd_solver/stencil_defines.h b/acc/mhd_solver/stencil_defines.h index 6e22bc9..26e7ffd 100644 --- a/acc/mhd_solver/stencil_defines.h +++ b/acc/mhd_solver/stencil_defines.h @@ -78,6 +78,7 @@ FUNC(AC_sink_pos_y),\ FUNC(AC_sink_pos_z),\ FUNC(AC_M_sink),\ + FUNC(AC_M_sink_init),\ FUNC(AC_M_sink_Msun),\ FUNC(AC_soft),\ FUNC(AC_accretion_range),\ diff --git a/src/standalone/config_loader.cc b/src/standalone/config_loader.cc index 95e4796..bee7050 100644 --- a/src/standalone/config_loader.cc +++ b/src/standalone/config_loader.cc @@ -150,6 +150,8 @@ update_config(AcMeshInfo* config) config->real_params[AC_M_sink] = config->real_params[AC_M_sink_Msun] * M_sun / config->real_params[AC_unit_mass]; + config->real_params[AC_M_sink_init] = config->real_params[AC_M_sink_Msun] * M_sun / + config->real_params[AC_unit_mass]; config->real_params[AC_G_const] = G_CONST_CGS / ((config->real_params[AC_unit_velocity] * config->real_params[AC_unit_velocity]) / diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 8159084..b9bf199 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -240,9 +240,16 @@ run_simulation(void) // acUpdate_sink_particle() will do the similar trick to the device. /* Step the simulation */ + AcReal accreted_mass = 0.0; for (int i = 1; i < max_steps; ++i) { const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); const AcReal dt = host_timestep(umax, mesh_info); + const AcReal sum_mass = acReduceScal(RTYPE_MAX, VTXBUF_ACCRETION); + accreted_mass = accreted_mass + sum_mass; + AcReal sink_mass = AC_M_sink_init + accreted_mass; + printf("sink mass is: %e\n", sink_mass); + printf("accreted mass is: %e\n", accreted_mass); + acLoadDeviceConstant(AC_M_sink, sink_mass); #if LFORCING const ForcingParams forcing_params = generateForcingParams(mesh_info);