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

Write now sink infor to time series. Also a bug cerrection for sink switch.

Browse Source
This commit is contained in:
Miikka Vaisala
2019-08-22 11:50:31 +08:00
parent 4414125a08
commit a81bc22fb6
4 changed files with 56 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
acc/mhd_solver/stencil_process.sps
View File

@@ -48,22 +48,27 @@ gradients(in VectorField uu)
#if LSINK
Vector
sink_gravity(int3 globalVertexIdx){
Vector force_gravity;
const Vector grid_pos = (Vector){(globalVertexIdx.x - nx_min) * dsx,
(globalVertexIdx.y - ny_min) * dsy,
(globalVertexIdx.z - nz_min) * dsz};
const Scalar sink_mass = DCONST_REAL(AC_M_sink);
const Vector sink_pos = (Vector){DCONST_REAL(AC_sink_pos_x),
DCONST_REAL(AC_sink_pos_y),
DCONST_REAL(AC_sink_pos_z)};
const Scalar distance = length(grid_pos - sink_pos);
const Scalar soft = DCONST_REAL(AC_soft);
const Scalar gravity_magnitude = (AC_G_const * sink_mass) / pow(((distance * distance) + soft*soft), 1.5);
const Vector direction = (Vector){(sink_pos.x - grid_pos.x) / distance,
(sink_pos.y - grid_pos.y) / distance,
(sink_pos.z - grid_pos.z) / distance};
force_gravity = gravity_magnitude * direction;
return force_gravity;
int accretion_switch = DCONST_INT(AC_switch_accretion);
if (accretion_switch == 1){
Vector force_gravity;
const Vector grid_pos = (Vector){(globalVertexIdx.x - nx_min) * dsx,
(globalVertexIdx.y - ny_min) * dsy,
(globalVertexIdx.z - nz_min) * dsz};
const Scalar sink_mass = DCONST_REAL(AC_M_sink);
const Vector sink_pos = (Vector){DCONST_REAL(AC_sink_pos_x),
DCONST_REAL(AC_sink_pos_y),
DCONST_REAL(AC_sink_pos_z)};
const Scalar distance = length(grid_pos - sink_pos);
const Scalar soft = DCONST_REAL(AC_soft);
const Scalar gravity_magnitude = (AC_G_const * sink_mass) / pow(((distance * distance) + soft*soft), 1.5);
const Vector direction = (Vector){(sink_pos.x - grid_pos.x) / distance,
(sink_pos.y - grid_pos.y) / distance,
(sink_pos.z - grid_pos.z) / distance};
force_gravity = gravity_magnitude * direction;
return force_gravity;
} else {
return (Vector){0.0, 0.0, 0.0};
}
}
#endif

17
analysis/python/astar/visual/lineplot.py
View File

@@ -32,7 +32,7 @@ def plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3):
plt.xlabel(xaxis)
plt.legend()
def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False, uuz=False, ss=False, acc=False):
def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False, uuz=False, ss=False, acc=False, sink=False):
if show_all:
lnrho=True
@@ -42,6 +42,7 @@ def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False,
uuz=True
ss=True
acc=True
sink=True
if lnrho:
plt.figure()
@@ -98,5 +99,19 @@ def plot_ts(ts, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False,
yaxis2 = 'acc_min'
yaxis3 = 'acc_max'
plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
plt.figure()
xaxis = 't_step'
yaxis1 = 'sink_mass'
yaxis2 = 'sink_mass'
yaxis3 = 'sink_mass'
plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
plt.figure()
xaxis = 't_step'
yaxis1 = 'accreted_mass'
yaxis2 = 'accreted_mass'
yaxis3 = 'accreted_mass'
plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3)
plt.show()

4
analysis/python/samples/readtest.py
View File

@@ -45,7 +45,7 @@ print("sys.argv", sys.argv)
#meshdir = "/tiara/home/mvaisala/astaroth-code/astaroth_2.0/build/"
#meshdir = "/tiara/ara/data/mvaisala/tmp/astaroth-code/astaroth_2.0/build/"
#meshdir = "/tiara/ara/data/mvaisala/asth_testbed_double/"
meshdir = "/scratch/data/mvaisala/forcingtest/"
meshdir = "/home/mvaisala/astaroth/build/"
if "xtopbound" in sys.argv:
for i in range(0, 171):
@@ -199,6 +199,6 @@ if 'sl' in sys.argv:
if 'ts' in sys.argv:
ts = ad.read.TimeSeries(fdir=meshdir)
vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, ss=1)
vis.lineplot.plot_ts(ts, show_all=True)

21
src/standalone/simulation.cc
View File

@@ -145,7 +145,8 @@ save_mesh(const AcMesh& save_mesh, const int step, const AcReal t_step)
// This function prints out the diagnostic values to std.out and also saves and
// appends an ascii file to contain all the result.
static inline void
print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* diag_file)
print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* diag_file,
const AcReal sink_mass, const AcReal accreted_mass)
{
AcReal buf_rms, buf_max, buf_min;
@@ -175,6 +176,10 @@ 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);
}
fprintf(diag_file, "\n");
}
@@ -213,11 +218,17 @@ run_simulation(void)
fprintf(diag_file, "%s_min %s_rms %s_max ", vtxbuf_names[i], vtxbuf_names[i],
vtxbuf_names[i]);
}
#if LSINK
fprintf(diag_file, "sink_mass accreted_mass ");
#endif
fprintf(diag_file, "\n");
write_mesh_info(&mesh_info);
print_diagnostics(0, AcReal(.0), t_step, diag_file);
#if LSINK
print_diagnostics(0, AcReal(.0), t_step, diag_file, mesh_info.real_params[AC_M_sink_init], 0.0);
#else
print_diagnostics(0, AcReal(.0), t_step, diag_file, -1.0, -1.0);
#endif
acBoundcondStep();
acStore(mesh);
@@ -265,6 +276,8 @@ run_simulation(void)
on_off_switch = 1;
}
acLoadDeviceConstant(AC_switch_accretion, on_off_switch);
#else
accreted_mass = -1.0; sink_mass = -1.0;
#endif
#if LFORCING
@@ -306,7 +319,7 @@ run_simulation(void)
timeseries.ts.
*/
print_diagnostics(i, dt, t_step, diag_file);
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);