diff --git a/examples/DiscPatch/HydroStatic/disc-patch-icc.yml b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml index 04d5b4cfd9ac45a218fd432db282258de2d151e4..9d1fe1b0948594f79c401e50308166313aa39924 100644 --- a/examples/DiscPatch/HydroStatic/disc-patch-icc.yml +++ b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml @@ -41,6 +41,6 @@ DiscPatchPotential: scale_height: 100. z_disc: 400. z_trunc: 300. - z_max: 380. + z_max: 350. timestep_mult: 0.03 growth_time: 5. diff --git a/examples/DiscPatch/HydroStatic/plotSolution.py b/examples/DiscPatch/HydroStatic/plotSolution.py index 845d1cd68c5118fdec2d7409519b7d3474352cee..e70f595c645eaadc282e511bbc0d2c510025cc87 100644 --- a/examples/DiscPatch/HydroStatic/plotSolution.py +++ b/examples/DiscPatch/HydroStatic/plotSolution.py @@ -35,6 +35,8 @@ import sys surface_density = 10. scale_height = 100. z_disc = 400. +z_trunc = 300. +z_max = 350. utherm = 20.2678457288 gamma = 5. / 3. @@ -78,22 +80,37 @@ for f in sorted(glob.glob("Disc-Patch_*.hdf5")): print "processing", f, "..." - zrange = np.linspace(0., 400., 1000) + zrange = np.linspace(0., 2. * z_disc, 1000) time, z, rho, P, v = get_data(f) fig, ax = pl.subplots(3, 1, sharex = True) ax[0].plot(z, rho, "r.") ax[0].plot(zrange, get_analytic_density(zrange), "k-") + ax[0].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5) + ax[0].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5) + ax[0].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5) + ax[0].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5) + ax[0].set_ylim(0., 1.2 * get_analytic_density(z_disc)) ax[0].set_ylabel("density") ax[1].plot(z, v, "r.") ax[1].plot(zrange, np.zeros(len(zrange)), "k-") + ax[1].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5) + ax[1].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5) + ax[1].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5) + ax[1].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5) + ax[1].set_ylim(-0.5, 10.) ax[1].set_ylabel("velocity norm") ax[2].plot(z, P, "r.") ax[2].plot(zrange, get_analytic_pressure(zrange), "k-") - ax[2].set_xlim(0., 400.) + ax[2].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5) + ax[2].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5) + ax[2].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5) + ax[2].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5) + ax[2].set_xlim(0., 2. * z_disc) + ax[2].set_ylim(0., 1.2 * get_analytic_pressure(z_disc)) ax[2].set_xlabel("z") ax[2].set_ylabel("pressure")