Coder Social home page Coder Social logo

swift_share_public's People

Contributors

jingyaodou avatar

Stargazers

avatar avatar

Watchers

avatar avatar avatar

Forkers

francyrad burgerz0

swift_share_public's Issues

Giant impact evolution in one plot

Dear @JingyaoDOU

A while ago I edited your sample script so that it could analyze the state of the particles after a giant impact, divided both by the color of the target particles and the impacting body. I must say that it works very well, I enclose an example here (sorry for the bad graphic).

Screenshot 2022-12-16 alle 01 29 26

The point is that it would be fabulous to make a video if the script could analyze a multitude of snapshots and save them in sequence. The point is, I don’t really know where to get my hands on it because there are a lot of definitions and I don’t really understand how I should load it into a new definition and then plot it.

Please, could you take a look? Where can i edit in way to get a sequence? In this example i deleted the old "swift load" because useless for the example, but maybe useful to plot a series of snapshot (?). Thank you very much.

Francesco

'
#!/usr/bin/env python

import h5py
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import swiftsimio as sw
import unyt
import woma
from woma.misc import utils
from copy import copy
R_earth = 6.371e6 # m
M_earth = 5.9724e24 # kg m^-3
G = 6.67408e-11 # m^3 kg^-1 s^-2

woma.load_eos_tables() # load the eos table for the calculation of entropy at later stage

def load_to_woma(snapshot):
# Load
data = sw.load(snapshot)
box_mid = 0.5 * data.metadata.boxsize[0].to(unyt.m)
data.gas.coordinates.convert_to_mks()
pos = np.array(data.gas.coordinates - box_mid)
data.gas.velocities.convert_to_mks()
vel = np.array(data.gas.velocities)
data.gas.smoothing_lengths.convert_to_mks()
h = np.array(data.gas.smoothing_lengths)
data.gas.masses.convert_to_mks()
m = np.array(data.gas.masses)
data.gas.densities.convert_to_mks()
rho = np.array(data.gas.densities)
data.gas.pressures.convert_to_mks()
p = np.array(data.gas.pressures)
data.gas.internal_energies.convert_to_mks()
u = np.array(data.gas.internal_energies)
matid = np.array(data.gas.material_ids)
# pid = np.array(data.gas.particle_ids)

pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
vel_centerM = np.sum(vel * m[:, np.newaxis], axis=0) / np.sum(m)

pos -= pos_centerM
vel -= vel_centerM

pos_noatmo = pos[matid != 200]

xy = np.hypot(pos_noatmo[:, 0], pos_noatmo[:, 1])
r = np.hypot(xy, pos_noatmo[:, 2])
r = np.sort(r)
R = np.mean(r[-1000:])

return pos, vel, h, m, rho, p, u, matid, R

Material IDs ( = type_id * type_factor + unit_id )

type_factor = 100
type_ANEOS = 4
type_Til = 1
type_HHe = 2
type_SESAME = 3
type_idg = 0
id_body = 200000000

Name and ID

Di_mat_id = {
"ANEOS_iron": type_ANEOS * type_factor + 1,
"ANEOS_iron_2": type_ANEOS * type_factor + 1 + id_body,
"ANEOS_alloy": type_ANEOS * type_factor + 2,
"ANEOS_alloy_2": type_ANEOS * type_factor + 2 + id_body,
"ANEOS_forsterite": type_ANEOS * type_factor,
"ANEOS_forsterite_2": type_ANEOS * type_factor + id_body,
"HM80_HHE": type_HHe * type_factor,
"SS08_water": type_SESAME * type_factor + 3,
"AQUA": type_SESAME * type_factor + 4,
"CMS19_HHe": type_SESAME * type_factor + 7,
"idg_HHe": type_idg * type_factor,
"Til_iron": type_Til * type_factor,
"Til_iron_2": type_Til * type_factor + id_body,
"Til_granite": type_Til * type_factor + 1,
"Til_granite_2": type_Til * type_factor + 1 + id_body,
}

Colour

Di_mat_colour = {
"ANEOS_iron": "black", # "tomato",
"ANEOS_alloy": "tomato",
"ANEOS_forsterite": "green", # "mediumseagreen",
"ANEOS_iron_2": "gold",
"ANEOS_alloy_2": "sandybrown",
"ANEOS_forsterite_2": "purple",
"SS08_water": "skyblue",
"AQUA": "skyblue",
"HM80_HHE": "aliceblue",
"CMS19_HHe": "lavenderblush",
"idg_HHe": "lavenderblush",
"Til_iron": "tomato",
"Til_iron_2": "sandybrown",
"Til_granite": "mediumseagreen",
"Til_granite_2": "pink",
}

marker size

Di_mat_size = {
"ANEOS_iron": 1,
"ANEOS_alloy": 1,
"ANEOS_forsterite": 1,
"ANEOS_alloy_2": 1,
"ANEOS_iron_2": 1,
"ANEOS_forsterite_2": 1,
"HM80_HHE": 1,
"SS08_water": 1,
"AQUA": 1,
"CMS19_HHe": 1,
"idg_HHe": 5,
"Til_iron": 1,
"Til_iron_2": 1,
"Til_granite": 1,
"Til_granite_2": 1,
}

Di_id_colour = {Di_mat_id[mat]: colour for mat, colour in Di_mat_colour.items()}

Di_id_size = {Di_mat_id[mat]: size * 3.5 for mat, size in Di_mat_size.items()}

density limits

rhomin = 5e-5
rhomax = 15.0

entropy limits

cmin = 1.5
cmax = 12.0

number of cells for grid

Ng = 601j
Ngz = 21j

region to grid/plot

xmax = 3
xmin = -xmax
ymin = -xmax
ymax = xmax
zmin = -3
zmax = 3

zcut = 6.0

def sw_plot(loc, snapshot_id=None, ax=None, npt=1000000, ax_lim=3.0, belowZ=True, figsz=10):
"""Select and load the particles to plot.

loc: snapshot path
snapshot_id: if loc is not the full path of a snapshot but a subparent path, then you need to provide the snapshot number you want to plot.
ax: the returned axis
npt: number of particles in the target
ax_lim: axis limit [-ax_lim,ax_lim]
belowZ: if only plot semisphere particles below the z=0 plane.
figsz: figure size
"""
# Load
if snapshot_id is not None:
    data = sw.load(loc + snapshot)
else:
    data = sw.load(loc)

box_mid = data.metadata.boxsize[0].to(unyt.Rearth)
pos = data.gas.coordinates.to(unyt.Rearth) - box_mid
id = data.gas.particle_ids
mat_id = data.gas.material_ids.value
m = data.gas.masses

# move to origin to the CoM
pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
pos -= pos_centerM

# Restrict to z < 0
if belowZ:
    sel = np.where(np.logical_and(pos[:, 2] < 0.1, pos[:, 2] > -600))[0]
    mat_id = mat_id[sel]
    id = id[sel]
    m = m[sel]
    pos = pos[sel]

# Sort in z order so higher particles are plotted on top
sort = np.argsort(pos[:, 2])
mat_id = mat_id[sort]
id = id[sort]
m = m[sort]
pos = pos[sort]

mat_id[npt <= id] += id_body

# Plot the particles, coloured by their material.
fig = plt.figure(figsize=(figsz, figsz))
if ax is not None:
    ax = ax
else:
    ax = plt.gca()

ax.set_aspect("equal", anchor="C")

colour = np.empty(len(pos), dtype=object)
for id_c, c in Di_id_colour.items():
    colour[mat_id == id_c] = c

size = np.empty(len(pos), dtype=float)
for id_s, s in Di_id_size.items():
    size[mat_id == id_s] = s

ax.scatter(
    pos[:, 0],
    pos[:, 1],
    s=3.5,
    c=colour,
    edgecolors="none",
    marker=".",
    alpha=1,
)

ax.set_ylabel(r"y Position ($R_\oplus$)", fontsize=22)
ax.set_xlabel(r"x Position ($R_\oplus$)", fontsize=22)

ax.yaxis.label.set_color("k")
ax.xaxis.label.set_color("k")
ax.tick_params(axis="x", colors="k", labelsize=14)
ax.tick_params(axis="y", colors="k", labelsize=14)

ax.set_xlim([-ax_lim, ax_lim])
ax.set_ylim([-ax_lim, ax_lim])

return ax

Ng = 1601j
Ngz = 201j

def sw_particles_plot(loc, ax_lim=2.0, Entropy=True, center=True, colourplot=True, npt=1000000, snapshot_id=None, figsz=10, ax=None):

"""plot the particles radial profile"""

data = sw.load(loc)
box_mid = 0.5 * data.metadata.boxsize[0].to(unyt.Rearth)
pos = data.gas.coordinates.to(unyt.Rearth) - box_mid
pos = np.array(pos)
data.gas.densities.convert_to_cgs()
rho_cgs = np.array(data.gas.densities)
data.gas.densities.convert_to_mks()
rho_mks = np.array(data.gas.densities)
data.gas.pressures.convert_to_cgs()
p = np.array(data.gas.pressures)
data.gas.internal_energies.convert_to_mks()
u = np.array(data.gas.internal_energies)
matid = data.gas.material_ids
data.gas.pressures.convert_to_mks()
p_mks = np.array(data.gas.pressures)
matid = np.array(matid)

id = data.gas.particle_ids
mat_id = data.gas.material_ids.value
m = data.gas.masses

colour = np.empty(len(matid), dtype=object)
for id_c, c in Di_id_colour.items():
    colour[matid == id_c] = c

entropy = np.zeros_like(p)
T = np.zeros_like(p)

mat_id[npt <= id] += id_body


if center:
    # data.gas.masses.convert_to_mks()
    m = np.array(data.gas.masses)
    pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
    # vel_centerM = np.sum(vel * m[:,np.newaxis], axis=0) / np.sum(m)

    pos -= pos_centerM
    # vel -= vel_centerM

if Entropy:
    sel = (matid != 200) & (matid != 0)
    entropy[sel] = woma.eos.eos.A1_s_u_rho(u[sel], rho_mks[sel], matid[sel])

T = woma.eos.eos.A1_T_u_rho(u, rho_mks, matid)

XY = np.hypot(pos[:, 0], pos[:, 1])
R = np.hypot(XY, pos[:, 2])





# Plot the particles, coloured by their material.
fig = plt.figure(figsize=(figsz, figsz))
if ax is not None:
    ax = ax
else:
    ax = plt.gca()

ax.set_aspect("equal", anchor="C")

colour = np.empty(len(pos), dtype=object)
for id_c, c in Di_id_colour.items():
    colour[mat_id == id_c] = c

size = np.empty(len(pos), dtype=float)
for id_s, s in Di_id_size.items():
    size[mat_id == id_s] = s

ax.scatter(
    pos[:, 0],
    pos[:, 1],
    s=3.5 * size,
    c=colour,
    edgecolors="none",
    marker=".",
    alpha=1,
)










print("Read %d particles!" % (len(R)))

fig, axs = plt.subplots(3, 2, figsize=(12, 15))
axs = axs.ravel()

axs[0] = sw_plot(loc, ax=axs[0])

if colourplot:
    axs[1].scatter(R, T, s=1, c=colour)
else:
    axs[1].scatter(R, T, s=1)
axs[1].tick_params(axis="x", direction="in")
axs[1].tick_params(axis="y", direction="in")
axs[1].set_title("Temperature profile", fontsize=22)
axs[1].set_xlabel(r"$R\oplus$")
axs[1].set_ylabel(r"T K")
axs[1].set_yscale("log")

if colourplot:
    axs[2].scatter(R, rho_cgs, s=1, c=colour)
else:
    axs[2].scatter(R, rho_cgs, s=1)
axs[2].set_title(r"Density $\rho$ profiles", fontsize=22)
axs[2].set_xlabel(r"$R\oplus$")
axs[2].set_ylabel(r"$\rho$ $g/cm^{3}$")
axs[2].set_yscale("log")
axs[2].tick_params(axis="x", direction="in")
axs[2].tick_params(axis="y", direction="in")
if colourplot:
    axs[3].scatter(R, p_mks / 1e9, s=1, c=colour)
else:
    axs[3].scatter(R, p_mks / 1e9, s=1)
axs[3].set_title("Pressure profile", fontsize=22)
axs[3].set_xlabel(r"$R\oplus$")
axs[3].set_ylabel("Gpa")
axs[3].set_yscale("log")
axs[3].tick_params(axis="x", direction="in")
axs[3].tick_params(axis="y", direction="in")

if colourplot:
    axs[4].scatter(R, entropy, s=1, c=colour)
else:
    axs[4].scatter(R, entropy, s=1)
axs[4].set_title("Entropy profile", fontsize=22)
axs[4].set_xlabel(r"$R\oplus$")
axs[4].set_ylabel("J/kg/K")
axs[4].tick_params(axis="x", direction="in")
axs[4].tick_params(axis="y", direction="in")

if colourplot:
    axs[5].scatter(R, u, s=1, c=colour)
else:
    axs[5].scatter(R, u, s=1)
axs[5].set_title("Internal energy profile", fontsize=22)
axs[5].set_xlabel(r"$R\oplus$")
axs[5].set_ylabel("J/kg")
axs[5].set_yscale("log")
axs[5].tick_params(axis="x", direction="in")
axs[5].tick_params(axis="y", direction="in")

fig.tight_layout()

plt.show()
plt.cla()
plt.clf()
plt.close()

sns.set_style("ticks")
sw_particles_plot("example.hdf5", colourplot=True)
'

Other scripts example

Dear Dr. Dao
Thank you so much for making the complete guide available, it will be very useful!
As you are a great expert on Python and the program, I ask you if it is possible to insert sample scripts where you can view other quantities of the simulation products instead of the material id, especially entropy and temperature (or internal energy, as your profile picture shows). Is it possible to have scripts that can see these quantities and save the snapshots of the simulation in sequence? In this way it would be possible to see the giant impacts also according to other physical quantities. Thank you very much for your availability!
Albert
internal_energy_example

initial conditions problems

Dear @JingyaoDOU

I wanted to unify 2 hdf5 files to make an impact initial conditions.

So, first of all, i created 2 planets and relaxed them, as you wrote.
I can read the planets both with swiftsimio and SPLASH, so, nothing wrong there.

These planets have fixed entropies enabled.

So, i created a script to unify them and it seems to work, there are no errors. However, when i try to start the simulation with swiftsim, it says that in the file there are no entropies data (i'm running swiftsim with fixed entropy as you suggest). Also, i cannot read the file with splash "impossible to read the header".

Did i do something wrong? Maybe I didn't load and appended all the data necessaries to make a collision? I attach you the file. Please, let me know
-Francesco

import h5py
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import swiftsimio as sw
import unyt
import woma
from woma.misc import utils
from copy import copy

R_earth = 6.371e6  # m
M_earth = 5.972e24  # kg m^-3

#The Gravitational constant, it may be useful if you want to add some constraints
G = 6.67408e-11  # m^3 kg^-1 s^-2

def load_to_woma(snapshot):
    # Load
    data = sw.load(snapshot)
    box_mid = 0.5 * data.metadata.boxsize[0].to(unyt.m)
    data.gas.coordinates.convert_to_mks()
    pos = np.array(data.gas.coordinates - box_mid)
    data.gas.velocities.convert_to_mks()
    vel = np.array(data.gas.velocities)
    data.gas.smoothing_lengths.convert_to_mks()
    h = np.array(data.gas.smoothing_lengths)
    data.gas.masses.convert_to_mks()
    m = np.array(data.gas.masses)
    data.gas.densities.convert_to_mks()
    rho = np.array(data.gas.densities)
    data.gas.pressures.convert_to_mks()
    p = np.array(data.gas.pressures)
    data.gas.internal_energies.convert_to_mks()
    u = np.array(data.gas.internal_energies)
    matid = np.array(data.gas.material_ids)
    # pid     = np.array(data.gas.particle_ids)

    pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
    vel_centerM = np.sum(vel * m[:, np.newaxis], axis=0) / np.sum(m)

    pos -= pos_centerM
    vel -= vel_centerM

    pos_noatmo = pos[matid != 200]

    xy = np.hypot(pos_noatmo[:, 0], pos_noatmo[:, 1])
    r = np.hypot(xy, pos_noatmo[:, 2])
    r = np.sort(r)
    R = np.mean(r[-1000:])

    return pos, vel, h, m, rho, p, u, matid, R
    
    
    
loc_tar = "Rufu_2017_target_1kk_1000.hdf5"
loc_imp = "Rufu_2017_impactor_1kk_0100.hdf5"

pos_tar, vel_tar, h_tar, m_tar, rho_tar, p_tar, u_tar, matid_tar, R_tar = load_to_woma(
    loc_tar)
    
pos_imp, vel_imp, h_imp, m_imp, rho_imp, p_imp, u_imp, matid_imp, R_imp = load_to_woma(
    loc_imp)

M_t = np.sum(m_tar)
M_i = np.sum(m_imp)
R_t = R_tar
R_i = R_imp


# Mutual escape speed
#v_esc = np.sqrt(2 * G * (M_target + M_impactor) / (R_target + R_impactor))

# Initial position and velocity of the target
A1_pos_t = np.array([0.0, 0.0, 0.0])
A1_vel_t = np.array([0.0, 0.0, 0.0])





A1_pos_i, A1_vel_i = woma.impact_pos_vel_b_v_c_r(
    #    b       = np.sin(45 * np.pi/180),
    b=0.7,
    v_c=20500,
    #t=3600,  # seconds. If you want time, you need to setup b_v_c_t.
    r=4*R_earth,
    R_t=R_t,
    R_i=R_i,
    M_t=M_t,
    M_i=M_i,
    #units_v_c="v_esc",
    #units_b="B"
    #b : float The impact parameter b=sin(B), or the impact angle B if units_b is "B".
)

A1_pos_com = (M_t * A1_pos_t + M_i * A1_pos_i) / (M_t + M_i)
A1_pos_t -= A1_pos_com
A1_pos_i -= A1_pos_com

# Centre of momentum
A1_vel_com = (M_t * A1_vel_t + M_i * A1_vel_i) / (M_t + M_i)
A1_vel_t -= A1_vel_com
A1_vel_i -= A1_vel_com

pos_tar += A1_pos_t
vel_tar[:] += A1_vel_t

pos_imp += A1_pos_i
vel_imp[:] += A1_vel_i
print("Number of particles in target %d" % (len(pos_tar)))


print(A1_pos_t / R_earth, "R_earth")
print(A1_vel_t, "m/s")
print(A1_pos_i / R_earth, "R_earth")
print(A1_vel_i, "m/s")



with h5py.File(
    "Rufu_est_momentum_1kk.hdf5", "w"
) as f:
    woma.save_particle_data(
        f,
        np.append(pos_tar, pos_imp, axis=0),
        np.append(vel_tar, vel_imp, axis=0),
        np.append(m_tar, m_imp),
        np.append(h_tar, h_imp),
        np.append(rho_tar, rho_imp),
        np.append(p_tar, p_imp),
        np.append(u_tar, u_imp),
        np.append(matid_tar, matid_imp),
        boxsize=100000 * R_earth,
        file_to_SI=woma.Conversions(M_earth, R_earth, 1),
    )

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.