diff --git a/src/standalone/model/host_memory.cc b/src/standalone/model/host_memory.cc index 1337a77..f6a54cb 100644 --- a/src/standalone/model/host_memory.cc +++ b/src/standalone/model/host_memory.cc @@ -225,9 +225,16 @@ simple_uniform_core(AcMesh* mesh) const double xorig = mesh->info.real_params[AC_xorig]; const double yorig = mesh->info.real_params[AC_yorig]; const double zorig = mesh->info.real_params[AC_zorig]; + + const double G_const = mesh->info.real_params[AC_G_const]; + const double M_sink_init = mesh->info.real_params[AC_M_sink_init]; + const double cs2_sound = mesh->info.real_params[AC_cs2_sound]; + + const double RR_inner_bound = mesh->info.real_params[AC_soft]/AcReal(2.0); + const double core_coeff = (exp(ampl_lnrho) * cs2_sound) / (double(4.0)*M_PI * G_const); + double xx, yy, zz, RR; - double delx, dely, delz; - double core_profile, core_coeff; + double core_profile; //TEMPORARY TEST INPUT PARAMETERS const double core_radius = DX*32.0; @@ -248,23 +255,15 @@ simple_uniform_core(AcMesh* mesh) yy = DY * double(j) - yorig; zz = DZ * double(k) - zorig; - delx = xx; - dely = yy; - delz = zz; - RR = sqrt(delx*delx + dely*dely + delz*delz); - - AcReal RR_inner_bound = mesh->info.real_params[AC_soft]/AcReal(2.0); - - core_coeff = (exp(ampl_lnrho) * mesh->info.real_params[AC_cs2_sound]) / - (double(4.0)*M_PI * mesh->info.real_params[AC_G_const]); + RR = sqrt(xx*xx + yy*yy + zz*zz); if (RR >= RR_inner_bound) { - abso_vel = vel_scale * sqrt(2.0 * mesh->info.real_params[AC_G_const] - * mesh->info.real_params[AC_M_sink_init] / RR); + abso_vel = vel_scale * sqrt(2.0 * G_const + * M_sink_init / RR); core_profile = pow(RR, -2.0); //double(1.0); } else { - abso_vel = vel_scale * sqrt(2.0 * mesh->info.real_params[AC_G_const] - * mesh->info.real_params[AC_M_sink_init] / RR_inner_bound); + abso_vel = vel_scale * sqrt(2.0 * AC_G_const + * AC_M_sink_init / RR_inner_bound); core_profile = pow(RR_inner_bound, -2.0); //double(1.0); } @@ -274,10 +273,10 @@ simple_uniform_core(AcMesh* mesh) } - mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(core_coeff*core_profile); - mesh->vertex_buffer[VTXBUF_UUX][idx] = -abso_vel * (yy / RR); - mesh->vertex_buffer[VTXBUF_UUY][idx] = abso_vel * (xx / RR); - mesh->vertex_buffer[VTXBUF_UUZ][idx] = double(0.0); + mesh->vertex_buffer[VTXBUF_LNRHO][idx] = AcReal(log(core_coeff*core_profile)); + mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(-abso_vel * (yy / RR)); + mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal( abso_vel * (xx / RR)); + mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(0.0); } } diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 012829c..efa6914 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -186,8 +186,8 @@ print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* di fprintf(diag_file, "%e %e %e ", double(buf_min), double(buf_rms), double(buf_max)); } - if ((sink_mass >= 0.0) || (accreted_mass >= 0.0)) { - fprintf(diag_file, "%e %e ", sink_mass, accreted_mass); + if ((sink_mass >= AcReal(0.0)) || (accreted_mass >= AcReal(0.0))) { + fprintf(diag_file, "%e %e ", double(sink_mass), double(accreted_mass)); } fprintf(diag_file, "\n"); @@ -208,8 +208,8 @@ run_simulation(void) AcMesh* mesh = acmesh_create(mesh_info); // TODO: This need to be possible to define in astaroth.conf - //acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, mesh); - acmesh_init_to(INIT_TYPE_SIMPLE_CORE, mesh); + acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, mesh); + //acmesh_init_to(INIT_TYPE_SIMPLE_CORE, mesh); //Initial condition for a collapse test #if LSINK vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh); @@ -285,6 +285,18 @@ run_simulation(void) on_off_switch = 1; } acLoadDeviceConstant(AC_switch_accretion, on_off_switch); + + //MV: Old TODOs to remind of eventual future directions. + //TODO_SINK acUpdate_sink_particle() + // 3. Velocity of the particle) + // TODO_SINK acAdvect_sink_particle() + // 1. Calculate the equation of motion for the sink particle. + // NOTE: Might require embedding with acIntegrate(dt). + // TODO_SINK acAccrete_sink_particle() + // 2. Transfer momentum into sink particle + // (OPTIONAL: Affection the motion of the particle) + // NOTE: Might require embedding with acIntegrate(dt). + // This is the hardest part. Please see Lee et al. ApJ 783 (2014) for reference. #else accreted_mass = -1.0; sink_mass = -1.0; #endif @@ -295,27 +307,7 @@ run_simulation(void) #endif - //TODO_SINK acUpdate_sink_particle() - // Update properties of the sing particle for acIntegrate(). Essentially: - // 1. Location of the particle - // 2. Mass of the particle - // 3. Velocity of the particle) - // These can be used for calculating he gravitational field. - // This is my first comment! by Jack - // This is my second comment! by Jack acIntegrate(dt); - // TODO_SINK acAdvect_sink_particle() - // THIS IS OPTIONAL. We will start from unmoving particle. - // 1. Calculate the equation of motion for the sink particle. - // NOTE: Might require embedding with acIntegrate(dt). - - // TODO_SINK acAccrete_sink_particle() - // Calculate accretion of the sink particle from the surrounding medium - // 1. Transfer density into sink particle mass - // 2. Transfer momentum into sink particle - // (OPTIONAL: Affection the motion of the particle) - // NOTE: Might require embedding with acIntegrate(dt). - // This is the hardest part. Please see Lee et al. ApJ 783 (2014) for reference. t_step += dt; @@ -329,9 +321,10 @@ run_simulation(void) */ print_diagnostics(i, dt, t_step, diag_file, sink_mass, accreted_mass); - printf("sink mass is: %.15e \n", sink_mass); - printf("accreted mass is: %.15e \n", accreted_mass); - +#if LSINK + printf("sink mass is: %.15e \n", double(sink_mass)); + printf("accreted mass is: %.15e \n", double(accreted_mass)); +#endif /* We would also might want an XY-average calculating funtion, which can be very useful when observing behaviour of turbulent