diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index e0c8211..d46c63a 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -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 diff --git a/analysis/python/astar/visual/lineplot.py b/analysis/python/astar/visual/lineplot.py index fe750ed..45d39df 100644 --- a/analysis/python/astar/visual/lineplot.py +++ b/analysis/python/astar/visual/lineplot.py @@ -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() diff --git a/analysis/python/samples/readtest.py b/analysis/python/samples/readtest.py index b2e4458..7f22858 100644 --- a/analysis/python/samples/readtest.py +++ b/analysis/python/samples/readtest.py @@ -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) diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index c16bdee..611c1c6 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -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);