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

Preparing isothermal collapse.

Browse Source
This commit is contained in:
Miikka Vaisala
2019-08-22 18:18:30 +08:00
parent a81bc22fb6
commit 1410e57866
8 changed files with 330 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
acc/mhd_solver/stencil_defines.h
View File

@@ -28,10 +28,10 @@
#define LDENSITY (1)
#define LHYDRO (1)
#define LMAGNETIC (1)
#define LENTROPY (1)
#define LENTROPY (0)
#define LTEMPERATURE (0)
#define LFORCING (1)
#define LUPWD (0)
#define LFORCING (0)
#define LUPWD (1)
#define LSINK (1)
#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter

8
acc/mhd_solver/stencil_process.sps
View File

@@ -238,7 +238,7 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala
}
#else
Vector
momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho) {
momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) {
Vector mom;
const Matrix S = stress_tensor(uu);
@@ -252,6 +252,10 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho) {
#if LSINK
+ sink_gravity(globalVertexIdx);
//Corresponding loss of momentum
- //(Scalar(1.0) / Scalar( (dsx*dsy*dsz) * exp(value(lnrho)))) * // Correction factor by unit mass
sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014)
;
#else
;
#endif
@@ -492,7 +496,7 @@ solve(Scalar dt) {
out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho, tt), dt);
out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt);
#else
out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho), dt);
out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho, dt), dt);
#endif
#if LFORCING

2
analysis/python/astar/data/read.py
View File

@@ -78,6 +78,8 @@ class Mesh:
if self.ok:
self.ss, timestamp, ok = read_bin('VTXBUF_ENTROPY', fdir, fnum, self.minfo)
self.accretion, timestamp, ok = read_bin('VTXBUF_ACCRETION', fdir, fnum, self.minfo)
#TODO Generalize is a dict. Do not hardcode!
uux, timestamp, ok = read_bin('VTXBUF_UUX', fdir, fnum, self.minfo)

247
analysis/python/jupyter/notebook_example.ipynb
View File

File diff suppressed because one or more lines are too long

79
src/standalone/model/host_memory.cc
View File

@@ -207,6 +207,81 @@ inflow_vedge(AcMesh* mesh)
}
}
// This is the initial condition type for the infalling vedge in the pseudodisk
// model.
void
simple_uniform_core(AcMesh* mesh)
{
const int mx = mesh->info.int_params[AC_mx];
const int my = mesh->info.int_params[AC_my];
const int mz = mesh->info.int_params[AC_mz];
//const int nx_min = mesh->info.int_params[AC_nx_min];
//const int nx_max = mesh->info.int_params[AC_nx_max];
//const int ny_min = mesh->info.int_params[AC_ny_min];
//const int ny_max = mesh->info.int_params[AC_ny_max];
//const int nz_min = mesh->info.int_params[AC_nz_min];
//const int nz_max = mesh->info.int_params[AC_nz_max];
const double DX = mesh->info.real_params[AC_dsx];
const double DY = mesh->info.real_params[AC_dsy];
const double DZ = mesh->info.real_params[AC_dsz];
const double ampl_lnrho = mesh->info.real_params[AC_ampl_lnrho];
//const double AMPL_UU = mesh->info.real_params[AC_ampl_uu];
// const double SQ2GM = mesh->info.real_params[AC_sq2GM_star];
// const double GM = mesh->info.real_params[AC_GM_star];
// const double M_star = mesh->info.real_params[AC_M_star];
// const double G_CONST = mesh->info.real_params[AC_G_CONST];
// const double unit_length = mesh->info.real_params[AC_unit_length];
// const double unit_density = mesh->info.real_params[AC_unit_density];
// const double unit_velocity = mesh->info.real_params[AC_unit_velocity];
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 trans = mesh->info.real_params[AC_trans];
double xx, yy, zz, RR;
double delx, dely, delz;
//double u_x, u_y, u_z;
double tanhRR;
//const double sink_pos_x = mesh->info.real_params[AC_sink_pos_x];
//const double sink_pos_y = mesh->info.real_params[AC_sink_pos_y];
//const double sink_pos_z = mesh->info.real_params[AC_sink_pos_z];
//TEMPORARY TEST INPUT PARAMETERS
const double core_radius = DX*32.0;
const double trans = DX*6.0;
for (int k = 0; k < mz; k++) {
for (int j = 0; j < my; j++) {
for (int i = 0; i < mx; i++) {
int idx = i + j*mx + k*mx*my;
xx = DX * double(i) - xorig;
yy = DY * double(j) - yorig;
zz = DZ * double(k) - zorig;
delx = xx;
dely = yy;
delz = zz;
RR = sqrt(delx*delx + dely*dely + delz*delz);
tanhRR = double(0.5)*tanh((core_radius-RR)/trans) + double(0.5)
+ double(0.1);
mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(exp(ampl_lnrho)*tanhRR);
mesh->vertex_buffer[VTXBUF_UUX][idx] = double(0.0);
mesh->vertex_buffer[VTXBUF_UUY][idx] = double(0.0);
mesh->vertex_buffer[VTXBUF_UUZ][idx] = double(0.0);
}
}
}
}
// This is the initial condition type for the infalling vedge in the pseudodisk
// model.
void
@@ -582,6 +657,10 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
}
}
break;
case INIT_TYPE_SIMPLE_CORE:
acmesh_clear(mesh);
simple_uniform_core(mesh);
break;
case INIT_TYPE_VEDGE:
acmesh_clear(mesh);
inflow_vedge_freefall(mesh);

1
src/standalone/model/host_memory.h
View File

@@ -34,6 +34,7 @@
FUNC(INIT_TYPE_XWAVE), \
FUNC(INIT_TYPE_GAUSSIAN_RADIAL_EXPL), \
FUNC(INIT_TYPE_ABC_FLOW) , \
FUNC(INIT_TYPE_SIMPLE_CORE), \
FUNC(INIT_TYPE_VEDGE), \
FUNC(INIT_TYPE_VEDGEX), \
FUNC(INIT_TYPE_RAYLEIGH_TAYLOR), \

5
src/standalone/renderer.cc
View File

@@ -292,7 +292,8 @@ renderer_quit(void)
return 0;
}
static int init_type = INIT_TYPE_GAUSSIAN_RADIAL_EXPL;
//static int init_type = INIT_TYPE_GAUSSIAN_RADIAL_EXPL;
static int init_type = INIT_TYPE_SIMPLE_CORE;
static bool
running(AcMesh* mesh)
@@ -396,7 +397,7 @@ run_renderer(void)
#if LSINK
const AcReal sum_mass = acReduceScal(RTYPE_MAX, VTXBUF_ACCRETION);
accreted_mass = accreted_mass + sum_mass;
AcReal sink_mass = AC_M_sink_init + accreted_mass;
AcReal sink_mass = mesh_info.real_params[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);

5
src/standalone/simulation.cc
View File

@@ -198,7 +198,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_GAUSSIAN_RADIAL_EXPL, mesh);
acmesh_init_to(INIT_TYPE_SIMPLE_CORE, mesh);
#if LSINK
vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh);
@@ -270,7 +271,7 @@ run_simulation(void)
vertex_buffer_set(VTXBUF_ACCRETION, 0.0, mesh);
int on_off_switch;
if (i < 1000) {
if (i < 1) {
on_off_switch = 0; //accretion is off till 1000 steps.
} else {
on_off_switch = 1;