Corrected type definition problems.
This commit is contained in:
@@ -225,9 +225,16 @@ simple_uniform_core(AcMesh* mesh)
|
|||||||
const double xorig = mesh->info.real_params[AC_xorig];
|
const double xorig = mesh->info.real_params[AC_xorig];
|
||||||
const double yorig = mesh->info.real_params[AC_yorig];
|
const double yorig = mesh->info.real_params[AC_yorig];
|
||||||
const double zorig = mesh->info.real_params[AC_zorig];
|
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 xx, yy, zz, RR;
|
||||||
double delx, dely, delz;
|
double core_profile;
|
||||||
double core_profile, core_coeff;
|
|
||||||
|
|
||||||
//TEMPORARY TEST INPUT PARAMETERS
|
//TEMPORARY TEST INPUT PARAMETERS
|
||||||
const double core_radius = DX*32.0;
|
const double core_radius = DX*32.0;
|
||||||
@@ -248,23 +255,15 @@ simple_uniform_core(AcMesh* mesh)
|
|||||||
yy = DY * double(j) - yorig;
|
yy = DY * double(j) - yorig;
|
||||||
zz = DZ * double(k) - zorig;
|
zz = DZ * double(k) - zorig;
|
||||||
|
|
||||||
delx = xx;
|
RR = sqrt(xx*xx + yy*yy + zz*zz);
|
||||||
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]);
|
|
||||||
|
|
||||||
if (RR >= RR_inner_bound) {
|
if (RR >= RR_inner_bound) {
|
||||||
abso_vel = vel_scale * sqrt(2.0 * mesh->info.real_params[AC_G_const]
|
abso_vel = vel_scale * sqrt(2.0 * G_const
|
||||||
* mesh->info.real_params[AC_M_sink_init] / RR);
|
* M_sink_init / RR);
|
||||||
core_profile = pow(RR, -2.0); //double(1.0);
|
core_profile = pow(RR, -2.0); //double(1.0);
|
||||||
} else {
|
} else {
|
||||||
abso_vel = vel_scale * sqrt(2.0 * mesh->info.real_params[AC_G_const]
|
abso_vel = vel_scale * sqrt(2.0 * AC_G_const
|
||||||
* mesh->info.real_params[AC_M_sink_init] / RR_inner_bound);
|
* AC_M_sink_init / RR_inner_bound);
|
||||||
core_profile = pow(RR_inner_bound, -2.0); //double(1.0);
|
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_LNRHO][idx] = AcReal(log(core_coeff*core_profile));
|
||||||
mesh->vertex_buffer[VTXBUF_UUX][idx] = -abso_vel * (yy / RR);
|
mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(-abso_vel * (yy / RR));
|
||||||
mesh->vertex_buffer[VTXBUF_UUY][idx] = abso_vel * (xx / RR);
|
mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal( abso_vel * (xx / RR));
|
||||||
mesh->vertex_buffer[VTXBUF_UUZ][idx] = double(0.0);
|
mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(0.0);
|
||||||
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@@ -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));
|
fprintf(diag_file, "%e %e %e ", double(buf_min), double(buf_rms), double(buf_max));
|
||||||
}
|
}
|
||||||
|
|
||||||
if ((sink_mass >= 0.0) || (accreted_mass >= 0.0)) {
|
if ((sink_mass >= AcReal(0.0)) || (accreted_mass >= AcReal(0.0))) {
|
||||||
fprintf(diag_file, "%e %e ", sink_mass, accreted_mass);
|
fprintf(diag_file, "%e %e ", double(sink_mass), double(accreted_mass));
|
||||||
}
|
}
|
||||||
|
|
||||||
fprintf(diag_file, "\n");
|
fprintf(diag_file, "\n");
|
||||||
@@ -208,8 +208,8 @@ run_simulation(void)
|
|||||||
|
|
||||||
AcMesh* mesh = acmesh_create(mesh_info);
|
AcMesh* mesh = acmesh_create(mesh_info);
|
||||||
// TODO: This need to be possible to define in astaroth.conf
|
// 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_GAUSSIAN_RADIAL_EXPL, mesh);
|
||||||
acmesh_init_to(INIT_TYPE_SIMPLE_CORE, mesh);
|
//acmesh_init_to(INIT_TYPE_SIMPLE_CORE, mesh); //Initial condition for a collapse test
|
||||||
|
|
||||||
#if LSINK
|
#if LSINK
|
||||||
vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh);
|
vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh);
|
||||||
@@ -285,6 +285,18 @@ run_simulation(void)
|
|||||||
on_off_switch = 1;
|
on_off_switch = 1;
|
||||||
}
|
}
|
||||||
acLoadDeviceConstant(AC_switch_accretion, on_off_switch);
|
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
|
#else
|
||||||
accreted_mass = -1.0; sink_mass = -1.0;
|
accreted_mass = -1.0; sink_mass = -1.0;
|
||||||
#endif
|
#endif
|
||||||
@@ -295,27 +307,7 @@ run_simulation(void)
|
|||||||
#endif
|
#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);
|
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;
|
t_step += dt;
|
||||||
|
|
||||||
@@ -329,9 +321,10 @@ run_simulation(void)
|
|||||||
*/
|
*/
|
||||||
|
|
||||||
print_diagnostics(i, dt, t_step, diag_file, sink_mass, accreted_mass);
|
print_diagnostics(i, dt, t_step, diag_file, sink_mass, accreted_mass);
|
||||||
printf("sink mass is: %.15e \n", sink_mass);
|
#if LSINK
|
||||||
printf("accreted mass is: %.15e \n", accreted_mass);
|
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,
|
We would also might want an XY-average calculating funtion,
|
||||||
which can be very useful when observing behaviour of turbulent
|
which can be very useful when observing behaviour of turbulent
|
||||||
|
Reference in New Issue
Block a user