Revising the initial condition into a self-similar profile.

This commit is contained in:
Miikka Vaisala
2019-09-03 18:42:14 +08:00
parent 65d69027be
commit 88a8198810
3 changed files with 21 additions and 11 deletions

View File

@@ -55,7 +55,7 @@ AC_sink_pos_y = 3.14
AC_sink_pos_z = 3.14
//AC_M_sink_Msun = 0.005
AC_M_sink_Msun = 1e-4 //1e-5
AC_M_sink_Msun = 1.0 // 1e-4 //1e-5
AC_soft = 0.25
AC_accretion_range = 0.25
@@ -71,5 +71,5 @@ AC_unit_length = 1.5e17
* Initial conditions
* =============================================================================
*/
AC_ampl_lnrho = 0.0
AC_ampl_uu = 0.5
AC_ampl_lnrho = 1.2
AC_ampl_uu = 1.0

View File

@@ -227,12 +227,12 @@ simple_uniform_core(AcMesh* mesh)
const double zorig = mesh->info.real_params[AC_zorig];
double xx, yy, zz, RR;
double delx, dely, delz;
double core_profile;
double core_profile, core_coeff;
//TEMPORARY TEST INPUT PARAMETERS
const double core_radius = DX*32.0;
const double trans = DX*12.0;
const double epsilon = DX*2.0;
//const double epsilon = DX*2.0;
const double vel_scale = mesh->info.real_params[AC_ampl_uu];
double abso_vel;
@@ -253,15 +253,19 @@ simple_uniform_core(AcMesh* mesh)
delz = zz;
RR = sqrt(delx*delx + dely*dely + delz*delz);
//core_profile = double(-0.5)*tanh((core_radius-RR)/trans) + double(0.5)
// + double(0.1);
core_profile = pow(RR+epsilon, -2.0); //double(1.0);
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) {
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 * mesh->info.real_params[AC_G_const]
* mesh->info.real_params[AC_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 * mesh->info.real_params[AC_G_const]
* mesh->info.real_params[AC_M_sink_init] / RR_inner_bound);
core_profile = pow(RR_inner_bound, -2.0); //double(1.0);
}
if (RR <= sqrt(DX*DX + DY*DY + DZ*DZ)) {
@@ -270,7 +274,7 @@ simple_uniform_core(AcMesh* mesh)
}
mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(exp(ampl_lnrho)*core_profile);
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);

View File

@@ -92,6 +92,12 @@ write_mesh_info(const AcMeshInfo* config)
fprintf(infotxt, "real AC_cs2_sound %e \n", (double)config->real_params[AC_cs2_sound]);
fprintf(infotxt, "real AC_cv_sound %e \n", (double)config->real_params[AC_cv_sound]);
//Physical units
fprintf(infotxt, "real AC_unit_density %e \n", (double)config->real_params[AC_unit_density]);
fprintf(infotxt, "real AC_unit_velocity %e \n", (double)config->real_params[AC_unit_velocity]);
fprintf(infotxt, "real AC_unit_mass %e \n", (double)config->real_params[AC_unit_mass ]);
fprintf(infotxt, "real AC_unit_length %e \n", (double)config->real_params[AC_unit_length ]);
//Here I'm still trying to copy the structure of the code above, and see if this will work for sink particle.
//I haven't fully undertand what these lines do but I'll read up on them soon. This is still yet experimental.
// Sink particle