|
| 1 | +import os.path |
| 2 | +from amuse.test.amusetest import get_path_to_results |
| 3 | +try: |
| 4 | + from matplotlib import pyplot |
| 5 | + HAS_MATPLOTLIB = True |
| 6 | + # from amuse.plot import plot, semilogy, xlabel, ylabel, loglog |
| 7 | +except ImportError: |
| 8 | + HAS_MATPLOTLIB = False |
| 9 | + |
| 10 | +from amuse.units import units |
| 11 | +from amuse.units import constants |
| 12 | +from amuse.units.generic_unit_converter import ConvertBetweenGenericAndSiUnits |
| 13 | +# from amuse.support.exceptions import AmuseException |
| 14 | +from amuse.community.mesa.interface import MESA |
| 15 | +# from amuse.community.gadget2.interface import Gadget2 |
| 16 | +from amuse.community.fi.interface import Fi |
| 17 | +from amuse.ext.star_to_sph import convert_stellar_model_to_SPH |
| 18 | + |
| 19 | +import numpy |
| 20 | + |
| 21 | +from amuse.datamodel import Particles |
| 22 | +from amuse.datamodel import ParticlesSuperset |
| 23 | +from amuse.datamodel import Grid |
| 24 | + |
| 25 | + |
| 26 | +def head_on_stellar_merger( |
| 27 | + masses=[0.3, 3.0] | units.MSun, |
| 28 | + star_age=310.0 | units.Myr, |
| 29 | + initial_separation=4.0 | units.RSun, |
| 30 | + angle=numpy.pi / 3, |
| 31 | + initial_speed=3000.0 | units.km / units.s, |
| 32 | + initial_speed_perpendicular=30.0 | units.km / units.s, |
| 33 | + number_of_sph_particles=50000, |
| 34 | + t_end=1.0e4 | units.s, |
| 35 | + sph_code=Fi, |
| 36 | + ): |
| 37 | + """ |
| 38 | + masses: Mass of the two stars |
| 39 | + star_age: Initial age of the stars |
| 40 | + number_of_sph_particles: Total number of particles of both stars, divided |
| 41 | + according to their masses |
| 42 | + t_end: (Physical, not computational) duration of the hydrodynamics |
| 43 | + simulation |
| 44 | + sph_code: Code to use for the hydrodynamics simulation |
| 45 | + """ |
| 46 | + |
| 47 | + # Convert some of the input parameters to string, for use in output file |
| 48 | + # names: |
| 49 | + n_string = "n" + ("%1.0e" % (number_of_sph_particles) |
| 50 | + ).replace("+0", "").replace("+", "") |
| 51 | + t_end_string = "t" + ("%1.0e" % (t_end.value_in(units.s)) |
| 52 | + ).replace("+0", "").replace("+", "") |
| 53 | + |
| 54 | + base_output_file_name = os.path.join( |
| 55 | + get_path_to_results(), "stellar_merger_"+n_string+"_"+t_end_string) |
| 56 | + |
| 57 | + stars = Particles(2) |
| 58 | + stars.mass = masses |
| 59 | + try: |
| 60 | + stellar_evolution = MESA() |
| 61 | + stellar_evolution.initialize_code() |
| 62 | + except: |
| 63 | + print("MESA was not built. Returning.") |
| 64 | + return |
| 65 | + stellar_evolution.commit_parameters() |
| 66 | + stellar_evolution.particles.add_particles(stars) |
| 67 | + stellar_evolution.commit_particles() |
| 68 | + print("Evolving stars with MESA...") |
| 69 | + stellar_evolution.evolve_model(star_age) |
| 70 | + |
| 71 | + number_of_sph_particles_1 = int( |
| 72 | + round( |
| 73 | + number_of_sph_particles |
| 74 | + * ( |
| 75 | + stellar_evolution.particles[0].mass |
| 76 | + / stellar_evolution.particles.mass.sum() |
| 77 | + ) |
| 78 | + ) |
| 79 | + ) |
| 80 | + number_of_sph_particles_2 = ( |
| 81 | + number_of_sph_particles - number_of_sph_particles_1 |
| 82 | + ) |
| 83 | + print("Creating initial conditions from a MESA stellar evolution model:") |
| 84 | + print( |
| 85 | + stellar_evolution.particles[0].mass, |
| 86 | + "star consisting of", |
| 87 | + number_of_sph_particles_1, |
| 88 | + "particles." |
| 89 | + ) |
| 90 | + sph_particles_1 = convert_stellar_model_to_SPH( |
| 91 | + stellar_evolution.particles[0], |
| 92 | + number_of_sph_particles_1, |
| 93 | + seed=12345 |
| 94 | + ).gas_particles |
| 95 | + print( |
| 96 | + stellar_evolution.particles[1].mass, |
| 97 | + "star consisting of", |
| 98 | + number_of_sph_particles_2, |
| 99 | + "particles." |
| 100 | + ) |
| 101 | + sph_particles_2 = convert_stellar_model_to_SPH( |
| 102 | + stellar_evolution.particles[1], |
| 103 | + number_of_sph_particles_2 |
| 104 | + ).gas_particles |
| 105 | + |
| 106 | + initial_separation += stellar_evolution.particles.radius.sum() |
| 107 | + sph_particles_2.x += numpy.cos(angle) * initial_separation |
| 108 | + sph_particles_2.y += numpy.sin(angle) * initial_separation |
| 109 | + sph_particles_1.vx += numpy.cos(angle) * initial_speed - \ |
| 110 | + numpy.sin(angle) * initial_speed_perpendicular |
| 111 | + sph_particles_1.vy += numpy.cos(angle) * initial_speed_perpendicular + \ |
| 112 | + numpy.sin(angle) * initial_speed |
| 113 | + view = [-0.5, 0.5, -0.5, 0.5] * \ |
| 114 | + (initial_separation + stellar_evolution.particles.radius.sum()) |
| 115 | + stellar_evolution.stop() |
| 116 | + |
| 117 | + all_sph_particles = ParticlesSuperset([sph_particles_1, sph_particles_2]) |
| 118 | + all_sph_particles.move_to_center() |
| 119 | + |
| 120 | + unit_converter = ConvertBetweenGenericAndSiUnits( |
| 121 | + 1.0 | units.RSun, constants.G, t_end) |
| 122 | + hydro_legacy_code = sph_code(unit_converter) |
| 123 | + n_steps = 100 |
| 124 | + hydro_legacy_code.parameters.n_smooth = 96 |
| 125 | + try: |
| 126 | + hydro_legacy_code.parameters.timestep = t_end / n_steps |
| 127 | + except Exception as exc: |
| 128 | + if "parameter is read-only" not in str(exc): |
| 129 | + raise |
| 130 | + hydro_legacy_code.gas_particles.add_particles(all_sph_particles) |
| 131 | + |
| 132 | + print("Evolving to t =", t_end, " (using", sph_code.__name__, "SPH code).") |
| 133 | + for time, i_step in [(i*t_end/n_steps, i) for i in range(1, n_steps+1)]: |
| 134 | + hydro_legacy_code.evolve_model(time) |
| 135 | + if not i_step % 4: |
| 136 | + hydro_plot( |
| 137 | + view, |
| 138 | + hydro_legacy_code, |
| 139 | + (300, 300), |
| 140 | + base_output_file_name + |
| 141 | + "_hydro_image{0:=03}.png".format(i_step) |
| 142 | + ) |
| 143 | + hydro_legacy_code.stop() |
| 144 | + print("All done!\n") |
| 145 | + |
| 146 | + |
| 147 | +def hydro_plot(view, hydro_code, image_size, figname): |
| 148 | + """ |
| 149 | + view: the (physical) region to plot [xmin, xmax, ymin, ymax] |
| 150 | + hydro_code: hydrodynamics code in which the gas to be plotted is defined |
| 151 | + image_size: size of the output image in pixels (x, y) |
| 152 | + """ |
| 153 | + shape = (image_size[0], image_size[1], 1) |
| 154 | + size = image_size[0] * image_size[1] |
| 155 | + axis_lengths = [0.0, 0.0, 0.0] | units.m |
| 156 | + axis_lengths[0] = view[1] - view[0] |
| 157 | + axis_lengths[1] = view[3] - view[2] |
| 158 | + grid = Grid.create(shape, axis_lengths) |
| 159 | + grid.x += view[0] |
| 160 | + grid.y += view[2] |
| 161 | + speed = grid.z.reshape(size) * (0 | 1/units.s) |
| 162 | + rho, rhovx, rhovy, rhovz, rhoe = \ |
| 163 | + hydro_code.get_hydro_state_at_point( |
| 164 | + grid.x.reshape(size), |
| 165 | + grid.y.reshape(size), |
| 166 | + grid.z.reshape(size), |
| 167 | + speed, speed, speed) |
| 168 | + |
| 169 | + min_v = 800.0 | units.km / units.s |
| 170 | + max_v = 3000.0 | units.km / units.s |
| 171 | + min_rho = 3.0e-4 | units.g / units.cm**3 |
| 172 | + max_rho = 0.1 | units.g / units.cm**3 |
| 173 | + min_E = 1.0e11 | units.J / units.kg |
| 174 | + max_E = 1.0e13 | units.J / units.kg |
| 175 | + |
| 176 | + v_sqr = (rhovx**2 + rhovy**2 + rhovz**2) / rho**2 |
| 177 | + E = rhoe / rho |
| 178 | + log_v = numpy.log((v_sqr / min_v**2)) / numpy.log((max_v**2 / min_v**2)) |
| 179 | + log_rho = numpy.log((rho / min_rho)) / numpy.log((max_rho / min_rho)) |
| 180 | + log_E = numpy.log((E / min_E)) / numpy.log((max_E / min_E)) |
| 181 | + |
| 182 | + red = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum( |
| 183 | + numpy.zeros_like(rho.number), log_rho)).reshape(shape) |
| 184 | + green = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum( |
| 185 | + numpy.zeros_like(rho.number), log_v)).reshape(shape) |
| 186 | + blue = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum( |
| 187 | + numpy.zeros_like(rho.number), log_E)).reshape(shape) |
| 188 | + alpha = numpy.minimum( |
| 189 | + numpy.ones_like(log_v), |
| 190 | + numpy.maximum( |
| 191 | + numpy.zeros_like(log_v), |
| 192 | + numpy.log((rho / (10*min_rho))) |
| 193 | + ) |
| 194 | + ).reshape(shape) |
| 195 | + |
| 196 | + rgba = numpy.concatenate((red, green, blue, alpha), axis=2) |
| 197 | + |
| 198 | + pyplot.figure(figsize=(image_size[0]/100.0, image_size[1]/100.0), dpi=100) |
| 199 | + im = pyplot.figimage(rgba, origin='lower') |
| 200 | + |
| 201 | + pyplot.savefig(figname, transparent=True, dpi = 100) |
| 202 | + print("\nHydroplot was saved to: ", figname) |
| 203 | + pyplot.close() |
| 204 | + |
| 205 | + |
| 206 | +if __name__ == "__main__": |
| 207 | + print("Running the simulation that formed the basis of the christmas card of Leiden Observatory of 2010.") |
| 208 | + print() |
| 209 | + print("Details:") |
| 210 | + print("The ornaments are the result of a smoothed particle simulation " \ |
| 211 | + "with 50000 equal mass particles of a 310 Myr star of 0.3 solar mass, " \ |
| 212 | + "which is ejected from a distance of 4 solar radii (left) " \ |
| 213 | + "with a velocity of 3000 km/s into a 3.0 solar mass star at an " \ |
| 214 | + "age of 310 Myr. The calculation was performed using the AMUSE " \ |
| 215 | + "(amusecode.org) software environment in which the stars were evolved " \ |
| 216 | + "using MESA to an age of 310 Myr before the encounter was performed " \ |
| 217 | + "using Fi. Each ornament, generated using pyplot, is a snapshot from the " \ |
| 218 | + "simulation, from top left to bottom right. The peak is created from a " \ |
| 219 | + "blend of all snapshots. The colors of the ornaments are, red: log " \ |
| 220 | + "of the density, green: log of the speed and for blue we used " \ |
| 221 | + "the log of the specific internal energy.") |
| 222 | + print() |
| 223 | + if HAS_MATPLOTLIB: |
| 224 | + head_on_stellar_merger() |
| 225 | + else: |
| 226 | + print("matplotlib is not installed. Install it in the site-packages folder of your Python installation. Returning.") |
0 commit comments