diff --git a/exemples/generate_vel_cube.py b/exemples/generate_vel_cube.py new file mode 100755 index 0000000000..af0d7c715a --- /dev/null +++ b/exemples/generate_vel_cube.py @@ -0,0 +1,133 @@ + +import numpy +import matplotlib.pyplot as plt + + + +def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): + from scipy import fftpack + + freqs = fftpack.fftfreq(res) + freq3d = numpy.array(numpy.meshgrid(freqs, freqs, freqs, indexing="ij")) + intfreq = numpy.around(freq3d * res) + kSqr = numpy.sum(numpy.abs(freq3d) ** 2, axis=0) + intkSqr = numpy.sum(numpy.abs(intfreq) ** 2, axis=0) + VK = [] + + # apply ~k^-2 exp(-k^2/kmax^2) filter to white noise to get x, y, and z components of velocity field + for i in range(3): + numpy.random.seed(seed + i) + rand_phase = fftpack.fftn( + numpy.random.normal(size=kSqr.shape) + ) # fourier transform of white noise + vk = rand_phase * (float(minmode) / res) ** 2 / (kSqr + 1e-300) + vk[intkSqr == 0] = 0.0 + vk[intkSqr < minmode**2] *= ( + intkSqr[intkSqr < minmode**2] ** 2 / minmode**4 + ) # smoother filter than mode-freezing; should give less "ringing" artifacts + vk *= numpy.exp(-intkSqr / maxmode**2) + + VK.append(vk) + VK = numpy.array(VK) + # bin = numpy.logspace(0,2.5,50) + # plt.hist(vk.flatten(),bins=bin) + # #plt.xlim(0,10**2.) + # plt.xscale("log") + # plt.yscale("log") + # plt.show() + # plt.imshow(vk[:,25,:].real) + # plt.show() + vk_new = numpy.zeros_like(VK) + + # do projection operator to get the correct mix of compressive and solenoidal + for i in range(3): + for j in range(3): + if i == j: + vk_new[i] += sol_weight * VK[j] + vk_new[i] += ( + (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j] + ) + vk_new[:, kSqr == 0] = 0.0 + VK = vk_new + + vel = numpy.array( + [fftpack.ifftn(vk).real for vk in VK] + ) # transform back to real space + vel -= numpy.average(vel, axis=(1, 2, 3))[:, numpy.newaxis, numpy.newaxis, numpy.newaxis] + vel = vel / numpy.sqrt(numpy.sum(vel**2, axis=0).mean()) # normalize so that RMS is 1 + return numpy.array(vel,dtype='f4') + + +# Global variables for velocity field +vx_global = None +vy_global = None +vz_global = None +domain_size_global = 1.0 + +def vel_field(pos): + """ + Interpolate velocity at position (x, y, z) using global velocity fields. + + Parameters: + ----------- + pos : tuple + (x, y, z) position at which to interpolate the velocity + + Returns: + -------- + tuple + (vx, vy, vz) velocity components at the given position(s) + """ + from scipy.interpolate import RegularGridInterpolator + + global vx_global, vy_global, vz_global, domain_size_global + + x, y, z = pos + res = vx_global.shape[0] + + # Create coordinate arrays for the grid + coords = numpy.linspace(0, domain_size_global, res) + + # Create interpolators for each velocity component + interp_vx = RegularGridInterpolator((coords, coords, coords), vx_global, + bounds_error=False, fill_value=0.0) + interp_vy = RegularGridInterpolator((coords, coords, coords), vy_global, + bounds_error=False, fill_value=0.0) + interp_vz = RegularGridInterpolator((coords, coords, coords), vz_global, + bounds_error=False, fill_value=0.0) + + points = numpy.column_stack([numpy.atleast_1d(x), + numpy.atleast_1d(y), + numpy.atleast_1d(z)]) + + vx_interp = interp_vx(points)[0] + vy_interp = interp_vy(points)[0] + vz_interp = interp_vz(points)[0] + + return vx_interp, vy_interp, vz_interp + + +if __name__ == "__main__": + + seed = 42 + + #avx,avy,avz = make_turb_field(res,power,seed) + vx,vy,vz = TurbField(128,2,64,0.7,seed) + + print(vx) + print(vy) + print(vz) + + # Set global velocity fields + vx_global = vx + vy_global = vy + vz_global = vz + domain_size_global = 1.0 + + # Example: Interpolate velocity at a specific position + test_pos = (0.5, 0.5, 0.5) + vel_x, vel_y, vel_z = vel_field(test_pos) + print(f"\nVelocity at {test_pos}:") + print(f" vx = {vel_x:.6f}") + print(f" vy = {vel_y:.6f}") + print(f" vz = {vel_z:.6f}") diff --git a/exemples/run_collapse.py b/exemples/run_collapse.py new file mode 100644 index 0000000000..cff417bf55 --- /dev/null +++ b/exemples/run_collapse.py @@ -0,0 +1,633 @@ +import glob +import json +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +import numpy + +import shamrock + +# Particle tracking is an experimental feature + +shamrock.enable_experimental_features() + +si = shamrock.UnitSystem() +sicte = shamrock.Constants(si) +codeu = shamrock.UnitSystem( + unit_time=sicte.year(), + unit_length=sicte.parsec(), + unit_mass=sicte.sol_mass(), +) +ucte = shamrock.Constants(codeu) + +m_s_codeu = codeu.get("m") * codeu.get("s", power=-1) +kg_m3_codeu = codeu.get("kg") * codeu.get("m", power=-3) + +print(f"m_s_codeu: {kg_m3_codeu}") + +cs = 190.0 # m/s +rho_g = kg_m3_codeu * 1.e-18 +initial_u = 1.0 +Npart = int(0.1e6) +ampl = 3 + +sim_radius = 2 + + +kb = ucte.kb() +print(f"kb: {kb}") +mu = 2.375 +mh = 1.00784 * ucte.dalton() +print(f"mu * mh * kb: {mu * mh * kb}") + +cs = cs * m_s_codeu # m/s +rho_c1 = 1.92e-13 * 1000 * kg_m3_codeu # g/cm^3 -> kg/m^3 +rho_c2 = 3.84e-8 * 1000 * kg_m3_codeu # g/cm^3 -> kg/m^3 +rho_c3 = 1.92e-3 * 1000 * kg_m3_codeu # g/cm^3 -> kg/m^3 + +sim_directory = f"collapse5_{Npart}" +import os + +os.makedirs(sim_directory, exist_ok=True) + + +tsound = sim_radius / cs +tff = shamrock.phys.free_fall_time(rho_g, codeu) + +P_init, _cs_init, T_init = shamrock.phys.eos.eos_Machida06( + cs=cs, rho=rho_g, rho_c1=rho_c1, rho_c2=rho_c2, rho_c3=rho_c3, mu=mu, mh=mh, kb=kb + ) + +print("---------------------------------------------------") +print(f"Npart : {Npart}") +print(f"sim_radius : {sim_radius} [pc]") +print(f"total mass : {rho_g * (2*sim_radius)**3:.3e} [solar mass]") +print(f"rho : {rho_g / kg_m3_codeu:.3e} [kg/m^3]") +print(f"P : {P_init / codeu.get("Pa"):.3e} [Pa]") +print(f"T : {T_init:.3e} [K]") +print(f"cs : {cs / m_s_codeu:.3e} [m/s]") +print(f"tsound (= R/cs) : {tsound:9.1f} [years]") +print(f"tff (= sqrt(3*pi/(32*G*rho))) : {tff:9.1f} [years]") +print(f"tff/tsound : {tff/tsound:.4f} (<1 = collapse)") +print("---------------------------------------------------") + + + + + + + + + + +class rho_xy_plot: + def __init__(self, model, center, ext_r, nx,ny, analysis_folder, analysis_prefix): + self.model = model + self.center = center + self.ext_r = ext_r + self.nx = nx + self.ny = ny + + self.analysis_prefix = os.path.join(analysis_folder, analysis_prefix) + "_" + self.plot_prefix = os.path.join(analysis_folder, "plot_" + analysis_prefix) + "_" + + self.npy_data_filename = self.analysis_prefix + "{:07}.npy" + self.json_data_filename = self.analysis_prefix + "{:07}.json" + self.json_glob = self.analysis_prefix + "*.json" + self.plot_filename = self.plot_prefix + "{:07}.png" + self.img_glob = self.analysis_prefix + "*.png" + + def compute_rho_xy(self): + + arr_rho_xy = model.render_cartesian_column_integ( + "rho", + "f64", + center=self.center, + delta_x=(self.ext_r * 2, 0, 0.0), + delta_y=(0.0, self.ext_r * 2, 0.0), + nx=self.nx, + ny=self.ny, + ) + + return arr_rho_xy + + def analysis_save(self, iplot): + arr_rho_xy = self.compute_rho_xy() + arr_rho_xy /= kg_m3_codeu + if shamrock.sys.world_rank() == 0: + cx,cy,cz = self.center + metadata = {"extent": [cx-self.ext_r, cx+self.ext_r, cy-self.ext_r, cy+self.ext_r], "time": self.model.get_time()} + np.save(self.npy_data_filename.format(iplot), arr_rho_xy) + + with open(self.json_data_filename.format(iplot), "w") as fp: + json.dump(metadata, fp) + + def load_analysis(self, iplot): + with open(self.json_data_filename.format(iplot), "r") as fp: + metadata = json.load(fp) + return np.load(self.npy_data_filename.format(iplot)), metadata + + def plot_rho_xy(self, iplot): + arr_rho_xy, metadata = self.load_analysis(iplot) + if shamrock.sys.world_rank() == 0: + + # Reset the figure using the same memory as the last one + plt.figure(num=1, clear=True, dpi=200) + + import copy + + my_cmap = matplotlib.colormaps["inferno"].copy() # copy the default cmap + my_cmap.set_bad(color="black") + + min_val = np.min(arr_rho_xy) + if min_val < 1e-20: + min_val = 1e-20 + + v_ext = 0.04 + res = plt.imshow( + arr_rho_xy, + cmap=my_cmap, + origin="lower", + extent=metadata["extent"], + norm="log", + vmin= min_val + ) + + plt.xlabel("x [pc]") + plt.ylabel("y [pc]") + plt.title("t = {:0.3f} [Year]".format(metadata["time"])) + + cbar = plt.colorbar(res, extend="both") + cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [kg.m$^{-2}$]") + + plt.savefig(self.plot_filename.format(iplot)) + plt.close() + + def get_list_dumps_id(self): + + list_files = glob.glob(self.json_glob) + list_files.sort() + list_json_files = [] + for f in list_files: + list_json_files.append(int(f.split("_")[-1].split(".")[0])) + return list_json_files + + def render_all(self): + if shamrock.sys.world_rank() == 0: + for iplot in self.get_list_dumps_id(): + print("Rendering rho xy plot for dump", iplot) + self.plot_rho_xy(iplot) + + + +class v_slice_xy_plot: + def __init__(self, model, center, ext_r, nx,ny, analysis_folder, analysis_prefix): + self.model = model + self.center = center + self.ext_r = ext_r + self.nx = nx + self.ny = ny + + self.analysis_prefix = os.path.join(analysis_folder, analysis_prefix) + "_" + self.plot_prefix = os.path.join(analysis_folder, "plot_" + analysis_prefix) + "_" + + self.npy_data_filename = self.analysis_prefix + "{:07}.npy" + self.json_data_filename = self.analysis_prefix + "{:07}.json" + self.json_glob = self.analysis_prefix + "*.json" + self.plot_filename = self.plot_prefix + "{:07}.png" + self.img_glob = self.analysis_prefix + "*.png" + + def compute_v_slice(self): + + arr_v_slice = model.render_cartesian_slice( + "vxyz", + "f64_3", + center=self.center, + delta_x=(self.ext_r * 2, 0, 0.0), + delta_y=(0.0, self.ext_r * 2, 0.0), + nx=self.nx, + ny=self.ny, + ) + + return arr_v_slice + + def analysis_save(self, iplot): + arr_v_slice = self.compute_v_slice() + arr_v_slice /= cs + if shamrock.sys.world_rank() == 0: + cx,cy,cz = self.center + metadata = {"extent": [cx-self.ext_r, cx+self.ext_r, cy-self.ext_r, cy+self.ext_r], "time": self.model.get_time()} + np.save(self.npy_data_filename.format(iplot), arr_v_slice) + + with open(self.json_data_filename.format(iplot), "w") as fp: + json.dump(metadata, fp) + + def load_analysis(self, iplot): + with open(self.json_data_filename.format(iplot), "r") as fp: + metadata = json.load(fp) + return np.load(self.npy_data_filename.format(iplot)), metadata + + def plot_v_slice(self, iplot): + arr_v_slice, metadata = self.load_analysis(iplot) + + v_norm = np.sqrt(arr_v_slice[:, :, 0] ** 2 + arr_v_slice[:, :, 1] ** 2 + arr_v_slice[:, :, 2] ** 2) + if shamrock.sys.world_rank() == 0: + + # Reset the figure using the same memory as the last one + plt.figure(num=1, clear=True, dpi=200) + + import copy + + my_cmap = matplotlib.colormaps["inferno"].copy() # copy the default cmap + my_cmap.set_bad(color="black") + + v_ext = 0.04 + res = plt.imshow( + v_norm, + cmap=my_cmap, + origin="lower", + extent=metadata["extent"] + ) + + plt.xlabel("x [pc]") + plt.ylabel("y [pc]") + plt.title("t = {:0.3f} [Year]".format(metadata["time"])) + + cbar = plt.colorbar(res, extend="both") + cbar.set_label(r"$||\mathbf{v}||/c_s$ ") + + plt.savefig(self.plot_filename.format(iplot)) + plt.close() + + def get_list_dumps_id(self): + + list_files = glob.glob(self.json_glob) + list_files.sort() + list_json_files = [] + for f in list_files: + list_json_files.append(int(f.split("_")[-1].split(".")[0])) + return list_json_files + + def render_all(self): + if shamrock.sys.world_rank() == 0: + for iplot in self.get_list_dumps_id(): + print("Rendering v slice plot for dump", iplot) + self.plot_v_slice(iplot) + + + + + +dump_prefix = sim_directory + "/dump_" +def get_dump_name(idump): + return dump_prefix + f"{idump:07}" + ".sham" + + +def get_last_dump(): + res = glob.glob(dump_prefix + "*.sham") + + num_max = -1 + + for f in res: + try: + dump_num = int(f[len(dump_prefix) : -5]) + if dump_num > num_max: + num_max = dump_num + except ValueError: + pass + + if num_max == -1: + return None + else: + return num_max + + +def purge_old_dumps(): + if shamrock.sys.world_rank() == 0: + res = glob.glob(dump_prefix + "*.sham") + res.sort() + + # The list of dumps to remove (keep the first and last 3 dumps) + to_remove = res[1:-3] + + for f in to_remove: + os.remove(f) + + + +bmin = (-sim_radius, -sim_radius, -sim_radius) +bmax = (sim_radius, sim_radius, sim_radius) + + +scheduler_split_val = int(2e7) +scheduler_merge_val = int(1) + + +N_target = Npart +xm, ym, zm = bmin +xM, yM, zM = bmax +vol_b = (xM - xm) * (yM - ym) * (zM - zm) + +if shamrock.sys.world_rank() == 0: + print("Npart", Npart) + print("scheduler_split_val", scheduler_split_val) + print("scheduler_merge_val", scheduler_merge_val) + print("N_target", N_target) + print("vol_b", vol_b) + +part_vol = vol_b / N_target + +# lattice volume +part_vol_lattice = 0.74 * part_vol + +dr = (part_vol_lattice / ((4.0 / 3.0) * 3.1416)) ** (1.0 / 3.0) + + +pmass = -1 + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4") + + +def is_in_sphere(pt): + x, y, z = pt + return (x**2 + y**2 + z**2) < sim_radius*sim_radius + +def config_model(): + cfg = model.gen_default_config() + # cfg.set_artif_viscosity_Constant(alpha_u = 1, alpha_AV = 1, beta_AV = 2) + # cfg.set_artif_viscosity_VaryingMM97(alpha_min = 0.1,alpha_max = 1,sigma_decay = 0.1, alpha_u = 1, beta_AV = 2) + cfg.set_artif_viscosity_VaryingCD10( + alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2 + ) + #cfg.set_boundary_periodic() + cfg.set_eos_machida06(rho_c1=rho_c1, rho_c2=rho_c2, rho_c3=rho_c3, cs=cs, mu=mu, mh=mh, kb=kb) + + cfg.set_self_gravity_fmm(order=5, opening_angle=0.5) + + cfg.set_units(codeu) + #cfg.print_status() + model.set_solver_config(cfg) + model.init_scheduler(scheduler_split_val, scheduler_merge_val) + + + bsim_min = (-sim_radius*2, -sim_radius*2, -sim_radius*2) + bsim_max = (sim_radius*2, sim_radius*2, sim_radius*2) + model.resize_simulation_box(bsim_min, bsim_max) + cfg.add_kill_sphere(center=(0, 0, 0), radius=sim_radius*2) # kill particles outside the simulation box + + setup = model.get_setup() + gen = setup.make_generator_lattice_hcp(dr, bmin, bmax) + + thesphere = setup.make_modifier_filter(parent=gen, filter=is_in_sphere) + + # On aurora /2 was correct to avoid out of memory + setup.apply_setup(thesphere, insert_step=int(scheduler_split_val / 2)) + + vol_b = (xM - xm) * (yM - ym) * (zM - zm) + totmass = rho_g * vol_b + pmass = model.total_mass_to_part_mass(totmass) + model.set_particle_mass(pmass) + + model.set_value_in_a_box("uint", "f64", initial_u, bmin, bmax) + + model.set_cfl_cour(0.1) + model.set_cfl_force(0.1) + + # if corrector is trigger CFL will become normal again after 1000 steps + model.set_cfl_mult_stiffness(5) + + +def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): + from scipy import fftpack + + freqs = fftpack.fftfreq(res) + freq3d = numpy.array(numpy.meshgrid(freqs, freqs, freqs, indexing="ij")) + intfreq = numpy.around(freq3d * res) + kSqr = numpy.sum(numpy.abs(freq3d) ** 2, axis=0) + intkSqr = numpy.sum(numpy.abs(intfreq) ** 2, axis=0) + VK = [] + + # apply ~k^-2 exp(-k^2/kmax^2) filter to white noise to get x, y, and z components of velocity field + for i in range(3): + numpy.random.seed(seed + i) + rand_phase = fftpack.fftn( + numpy.random.normal(size=kSqr.shape) + ) # fourier transform of white noise + vk = rand_phase * (float(minmode) / res) ** 2 / (kSqr + 1e-300) + vk[intkSqr == 0] = 0.0 + vk[intkSqr < minmode**2] *= ( + intkSqr[intkSqr < minmode**2] ** 2 / minmode**4 + ) # smoother filter than mode-freezing; should give less "ringing" artifacts + vk *= numpy.exp(-intkSqr / maxmode**2) + + VK.append(vk) + VK = numpy.array(VK) + # bin = numpy.logspace(0,2.5,50) + # plt.hist(vk.flatten(),bins=bin) + # #plt.xlim(0,10**2.) + # plt.xscale("log") + # plt.yscale("log") + # plt.show() + # plt.imshow(vk[:,25,:].real) + # plt.show() + vk_new = numpy.zeros_like(VK) + + # do projection operator to get the correct mix of compressive and solenoidal + for i in range(3): + for j in range(3): + if i == j: + vk_new[i] += sol_weight * VK[j] + vk_new[i] += ( + (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j] + ) + vk_new[:, kSqr == 0] = 0.0 + VK = vk_new + + vel = numpy.array( + [fftpack.ifftn(vk).real for vk in VK] + ) # transform back to real space + vel -= numpy.average(vel, axis=(1, 2, 3))[:, numpy.newaxis, numpy.newaxis, numpy.newaxis] + vel = vel / numpy.sqrt(numpy.sum(vel**2, axis=0).mean()) # normalize so that RMS is 1 + return numpy.array(vel,dtype='f4') + + + + +seed = 42 + +#avx,avy,avz = make_turb_field(res,power,seed) +vx,vy,vz = TurbField(128,2,64,0.7,seed) + +print(f"vx shape: {vx.shape}, dtype: {vx.dtype}") +print(f"vy shape: {vy.shape}, dtype: {vy.dtype}") +print(f"vz shape: {vz.shape}, dtype: {vz.dtype}") +print(f"vx min/max: {vx.min():.6f}/{vx.max():.6f}") +print(f"vy min/max: {vy.min():.6f}/{vy.max():.6f}") +print(f"vz min/max: {vz.min():.6f}/{vz.max():.6f}") + +# Set global velocity fields +vx_global = vx +vy_global = vy +vz_global = vz +domain_size_global = 1.0 +sim_bmin = bmin # (-0.5, -0.5, -0.5) +sim_bmax = bmax # (0.5, 0.5, 0.5) + +res = vx_global.shape[0] + +# Create coordinate arrays for the grid +coords = numpy.linspace(0, domain_size_global, res) + +# Create interpolators for each velocity component +from scipy.interpolate import RegularGridInterpolator +interp_vx = RegularGridInterpolator((coords, coords, coords), vx_global, + bounds_error=False, fill_value=0.0) +interp_vy = RegularGridInterpolator((coords, coords, coords), vy_global, + bounds_error=False, fill_value=0.0) +interp_vz = RegularGridInterpolator((coords, coords, coords), vz_global, + bounds_error=False, fill_value=0.0) + + +def vel_field(pos): + """ + Interpolate velocity at position (x, y, z) using global velocity fields. + + Parameters: + ----------- + pos : tuple + (x, y, z) position in simulation coordinates + + Returns: + -------- + tuple + (vx, vy, vz) velocity components at the given position(s) + """ + + + global interp_vx, interp_vy, interp_vz, domain_size_global, sim_bmin, sim_bmax, ampl, cs + + x, y, z = pos + + # Transform from simulation coordinates to velocity field coordinates [0, 1] + x_vf = (x - sim_bmin[0]) / (sim_bmax[0] - sim_bmin[0]) + y_vf = (y - sim_bmin[1]) / (sim_bmax[1] - sim_bmin[1]) + z_vf = (z - sim_bmin[2]) / (sim_bmax[2] - sim_bmin[2]) + + + + points = numpy.column_stack([numpy.atleast_1d(x_vf), + numpy.atleast_1d(y_vf), + numpy.atleast_1d(z_vf)]) + + vx_interp = interp_vx(points)[0] + vy_interp = interp_vy(points)[0] + vz_interp = interp_vz(points)[0] + + #print(f"vx_interp = {vx_interp}, vy_interp = {vy_interp}, vz_interp = {vz_interp}") + + return vx_interp*ampl *cs, vy_interp*ampl *cs, vz_interp*ampl *cs + +# Example: Interpolate velocity at a specific position (in simulation coordinates) +test_pos = (0.0, 0.0, 0.0) # Center of simulation domain +vel_x, vel_y, vel_z = vel_field(test_pos) +print(f"\nVelocity at {test_pos} (sim coords):") +print(f" vx = {vel_x:.6f}") +print(f" vy = {vel_y:.6f}") +print(f" vz = {vel_z:.6f}") + +# Test at edge of domain +test_pos2 = (0.25, -0.25, 0.1) +vel_x2, vel_y2, vel_z2 = vel_field(test_pos2) +print(f"\nVelocity at {test_pos2} (sim coords):") +print(f" vx = {vel_x2:.6f}") +print(f" vy = {vel_y2:.6f}") +print(f" vz = {vel_z2:.6f}") + +def setup_particles(): + + model.set_field_value_lambda_f64_3("vxyz", vel_field) + + model.timestep() + + +rho_plot_large = rho_xy_plot(model, (0.0, 0.0, 0.0), sim_radius *1.2, 2048, 2048, sim_directory, "rho_xy") +rho_plot_mid = rho_xy_plot(model, (0.0, 0.0, 0.0), sim_radius /1.2, 2048, 2048, sim_directory, "rho_mid") +rho_plot_zoom = rho_xy_plot(model, (0.0, 0.0, 0.0), 0.25, 2048, 2048, sim_directory, "rho_zoom") + +v_slice_plot = v_slice_xy_plot(model, (0.0, 0.0, 0.0), sim_radius *1.2, 2048, 2048, sim_directory, "v_slice") + +def analysis(i): + + rho_plot_large.analysis_save(i) + rho_plot_mid.analysis_save(i) + v_slice_plot.analysis_save(i) + + dat = ctx.collect_data() + + #get index with max rho + max_rho_index = np.argmin(dat["hpart"]) + hpart_min = dat["hpart"][max_rho_index] + max_rho = model.get_particle_mass() * (model.get_hfact() / hpart_min) ** 3 + + max_rho_x = dat["xyz"][max_rho_index, 0] + max_rho_y = dat["xyz"][max_rho_index, 1] + max_rho_z = dat["xyz"][max_rho_index, 2] + print(f"max rho: {max_rho/kg_m3_codeu:.3e} at ({max_rho_x:.3e}, {max_rho_y:.3e}, {max_rho_z:.3e})") + + #render around max rho + rho_plot_zoom.center = (max_rho_x, max_rho_y, max_rho_z) + + rho_plot_zoom.analysis_save(i) + + +idump_last_dump = get_last_dump() + +if shamrock.sys.world_rank() == 0: + print("Last dump:", idump_last_dump) + +# %% +# Load the last dump if it exists, setup otherwise + +run_simulation = "1" +i = 0 +next_t = 0 + +if run_simulation == "1" and idump_last_dump is not None: + model.load_from_dump(get_dump_name(idump_last_dump)) + i = idump_last_dump + 1 + next_t = model.get_time() + tff * 0.1 +elif run_simulation == "1": + config_model() + setup_particles() + + +rho_plot_large.render_all() +rho_plot_mid.render_all() +rho_plot_zoom.render_all() +v_slice_plot.render_all() + +while i < 10000: + + model.evolve_until(next_t, niter_max=250) + + next_t = model.get_time() + tff * 1e-2 + + analysis(i) + model.dump(get_dump_name(i)) + purge_old_dumps() + + if i % 2 == 0: + rho_plot_large.render_all() + rho_plot_mid.render_all() + rho_plot_zoom.render_all() + v_slice_plot.render_all() + i += 1 + + +rho_plot_large.render_all() +rho_plot_mid.render_all() +rho_plot_zoom.render_all() +v_slice_plot.render_all() \ No newline at end of file diff --git a/src/shammath/include/shammath/Perlin.h b/src/shammath/include/shammath/Perlin.h new file mode 100644 index 0000000000..a0cc7577b6 --- /dev/null +++ b/src/shammath/include/shammath/Perlin.h @@ -0,0 +1,370 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2025 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +#pragma once + +/** + * @file paving_function.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/sycl.hpp" + +namespace shammath { + + class PerlinNoise { + + static constexpr int grad3[][3] + = {{1, 1, 0}, + {-1, 1, 0}, + {1, -1, 0}, + {-1, -1, 0}, + {1, 0, 1}, + {-1, 0, 1}, + {1, 0, -1}, + {-1, 0, -1}, + {0, 1, 1}, + {0, -1, 1}, + {0, 1, -1}, + {0, -1, -1}}; + + inline double dot(const int *g, double x, double y) { return g[0] * x + g[1] * y; } + + inline double dot(const int *g, double x, double y, double z) { + return g[0] * x + g[1] * y + g[2] * z; + } + + inline double dot(const int *g, double x, double y, double z, double w) { + return g[0] * x + g[1] * y + g[2] * z + g[3] * w; + } + static constexpr int grad4[][4] + = {{0, 1, 1, 1}, {0, 1, 1, -1}, {0, 1, -1, 1}, {0, 1, -1, -1}, {0, -1, 1, 1}, + {0, -1, 1, -1}, {0, -1, -1, 1}, {0, -1, -1, -1}, {1, 0, 1, 1}, {1, 0, 1, -1}, + {1, 0, -1, 1}, {1, 0, -1, -1}, {-1, 0, 1, 1}, {-1, 0, 1, -1}, {-1, 0, -1, 1}, + {-1, 0, -1, -1}, {1, 1, 0, 1}, {1, 1, 0, -1}, {1, -1, 0, 1}, {1, -1, 0, -1}, + {-1, 1, 0, 1}, {-1, 1, 0, -1}, {-1, -1, 0, 1}, {-1, -1, 0, -1}, {1, 1, 1, 0}, + {1, 1, -1, 0}, {1, -1, 1, 0}, {1, -1, -1, 0}, {-1, 1, 1, 0}, {-1, 1, -1, 0}, + {-1, -1, 1, 0}, {-1, -1, -1, 0}}; + + static constexpr int p[] + = {151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, + 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, + 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, + 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, + 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, + 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, + 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, + 200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, + 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212, + 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213, + 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9, + 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, + 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, + 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, + 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, 205, 93, + 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180}; + + inline int perm(int i) { return p[i & 255]; } + + static constexpr int simplex[][4] + = {{0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, + {0, 0, 0, 0}, {1, 2, 3, 0}, {0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1}, + {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, + {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, + {1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, + {2, 3, 0, 1}, {2, 3, 1, 0}, {1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, + {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, + {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, + {2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {3, 0, 1, 2}, {3, 0, 2, 1}, + {0, 0, 0, 0}, {3, 1, 2, 0}, {2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, + {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}}; + + inline int fastfloor(double x) { return x > 0 ? (int) x : (int) x - 1; } + + public: + inline double noise(double xin, double yin) { + double n0, n1, n2; + + double F2 = 0.5 * (sycl::sqrt(3.0) - 1.0); + double s = (xin + yin) * F2; + + int i = fastfloor(xin + s); + int j = fastfloor(yin + s); + double G2 = (3.0 - sycl::sqrt(3.0)) / 6.0; + double t = (i + j) * G2; + double X0 = i - t; + + double Y0 = j - t; + double x0 = xin - X0; + + double y0 = yin - Y0; + int i1, j1; + if (x0 > y0) { + i1 = 1; + j1 = 0; + } else { + i1 = 0; + j1 = 1; + } + double x1 = x0 - i1 + G2; + + double y1 = y0 - j1 + G2; + double x2 = x0 - 1.0 + 2.0 * G2; + double y2 = y0 - 1.0 + 2.0 * G2; + int ii = i & 255; + int jj = j & 255; + int gi0 = perm(ii + perm(jj)) % 12; + int gi1 = perm(ii + i1 + perm(jj + j1)) % 12; + int gi2 = perm(ii + 1 + perm(jj + 1)) % 12; + double t0 = 0.5 - x0 * x0 - y0 * y0; + if (t0 < 0) + n0 = 0.0; + else { + t0 *= t0; + n0 = t0 * t0 * dot(grad3[gi0], x0, y0); + } + double t1 = 0.5 - x1 * x1 - y1 * y1; + if (t1 < 0) + n1 = 0.0; + else { + t1 *= t1; + n1 = t1 * t1 * dot(grad3[gi1], x1, y1); + } + + double t2 = 0.5 - x2 * x2 - y2 * y2; + if (t2 < 0) + n2 = 0.0; + else { + t2 *= t2; + n2 = t2 * t2 * dot(grad3[gi2], x2, y2); + } + + return 70.0 * (n0 + n1 + n2); + } + + inline double noise_3d(double xin, double yin, double zin) { + double n0, n1, n2, n3; + double F3 = 1.0 / 3.0; + double s = (xin + yin + zin) * F3; + int i = fastfloor(xin + s); + int j = fastfloor(yin + s); + int k = fastfloor(zin + s); + double G3 = 1.0 / 6.0; + double t = (i + j + k) * G3; + double X0 = i - t; + double Y0 = j - t; + double Z0 = k - t; + double x0 = xin - X0; + double y0 = yin - Y0; + double z0 = zin - Z0; + int i1, j1, k1; + int i2, j2, k2; + if (x0 >= y0) { + if (y0 >= z0) { + i1 = 1; + j1 = 0; + k1 = 0; + i2 = 1; + j2 = 1; + k2 = 0; + } else if (x0 >= z0) { + i1 = 1; + j1 = 0; + k1 = 0; + i2 = 1; + j2 = 0; + k2 = 1; + } else { + i1 = 0; + j1 = 0; + k1 = 1; + i2 = 1; + j2 = 0; + k2 = 1; + } + } else { + if (y0 < z0) { + i1 = 0; + j1 = 0; + k1 = 1; + i2 = 0; + j2 = 1; + k2 = 1; + } else if (x0 < z0) { + i1 = 0; + j1 = 1; + k1 = 0; + i2 = 0; + j2 = 1; + k2 = 1; + } else { + i1 = 0; + j1 = 1; + k1 = 0; + i2 = 1; + j2 = 1; + k2 = 0; + } + } + + double x1 = x0 - i1 + G3; + double y1 = y0 - j1 + G3; + double z1 = z0 - k1 + G3; + double x2 = x0 - i2 + 2.0 * G3; + double y2 = y0 - j2 + 2.0 * G3; + double z2 = z0 - k2 + 2.0 * G3; + double x3 = x0 - 1.0 + 3.0 * G3; + double y3 = y0 - 1.0 + 3.0 * G3; + double z3 = z0 - 1.0 + 3.0 * G3; + int ii = i & 255; + int jj = j & 255; + int kk = k & 255; + int gi0 = perm(ii + perm(jj + perm(kk))) % 12; + int gi1 = perm(ii + i1 + perm(jj + j1 + perm(kk + k1))) % 12; + int gi2 = perm(ii + i2 + perm(jj + j2 + perm(kk + k2))) % 12; + int gi3 = perm(ii + 1 + perm(jj + 1 + perm(kk + 1))) % 12; + double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0; + if (t0 < 0) + n0 = 0.0; + else { + t0 *= t0; + n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0); + } + double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1; + if (t1 < 0) + n1 = 0.0; + else { + t1 *= t1; + n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1); + } + double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2; + if (t2 < 0) + n2 = 0.0; + else { + t2 *= t2; + n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2); + } + double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3; + if (t3 < 0) + n3 = 0.0; + else { + t3 *= t3; + n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3); + } + return 32.0 * (n0 + n1 + n2 + n3); + } + + inline double noise(double x, double y, double z, double w) { + double F4 = (sycl::sqrt(5.0) - 1.0) / 4.0; + double G4 = (5.0 - sycl::sqrt(5.0)) / 20.0; + double n0, n1, n2, n3, n4; + double s = (x + y + z + w) * F4; + int i = fastfloor(x + s); + int j = fastfloor(y + s); + int k = fastfloor(z + s); + int l = fastfloor(w + s); + double t = (i + j + k + l) * G4; + double X0 = i - t; + double Y0 = j - t; + double Z0 = k - t; + double W0 = l - t; + double x0 = x - X0; + double y0 = y - Y0; + double z0 = z - Z0; + double w0 = w - W0; + int c1 = (x0 > y0) ? 32 : 0; + int c2 = (x0 > z0) ? 16 : 0; + int c3 = (y0 > z0) ? 8 : 0; + int c4 = (x0 > w0) ? 4 : 0; + int c5 = (y0 > w0) ? 2 : 0; + int c6 = (z0 > w0) ? 1 : 0; + int c = c1 + c2 + c3 + c4 + c5 + c6; + int i1, j1, k1, l1; + int i2, j2, k2, l2; + int i3, j3, k3, l3; + i1 = simplex[c][0] >= 3 ? 1 : 0; + j1 = simplex[c][1] >= 3 ? 1 : 0; + k1 = simplex[c][2] >= 3 ? 1 : 0; + l1 = simplex[c][3] >= 3 ? 1 : 0; + i2 = simplex[c][0] >= 2 ? 1 : 0; + j2 = simplex[c][1] >= 2 ? 1 : 0; + + k2 = simplex[c][2] >= 2 ? 1 : 0; + l2 = simplex[c][3] >= 2 ? 1 : 0; + i3 = simplex[c][0] >= 1 ? 1 : 0; + j3 = simplex[c][1] >= 1 ? 1 : 0; + k3 = simplex[c][2] >= 1 ? 1 : 0; + l3 = simplex[c][3] >= 1 ? 1 : 0; + double x1 = x0 - i1 + G4; + double y1 = y0 - j1 + G4; + double z1 = z0 - k1 + G4; + double w1 = w0 - l1 + G4; + double x2 = x0 - i2 + 2.0 * G4; + double y2 = y0 - j2 + 2.0 * G4; + double z2 = z0 - k2 + 2.0 * G4; + double w2 = w0 - l2 + 2.0 * G4; + double x3 = x0 - i3 + 3.0 * G4; + double y3 = y0 - j3 + 3.0 * G4; + double z3 = z0 - k3 + 3.0 * G4; + double w3 = w0 - l3 + 3.0 * G4; + double x4 = x0 - 1.0 + 4.0 * G4; + double y4 = y0 - 1.0 + 4.0 * G4; + double z4 = z0 - 1.0 + 4.0 * G4; + double w4 = w0 - 1.0 + 4.0 * G4; + int ii = i & 255; + int jj = j & 255; + int kk = k & 255; + int ll = l & 255; + int gi0 = perm(ii + perm(jj + perm(kk + perm(ll)))) % 32; + int gi1 = perm(ii + i1 + perm(jj + j1 + perm(kk + k1 + perm(ll + l1)))) % 32; + int gi2 = perm(ii + i2 + perm(jj + j2 + perm(kk + k2 + perm(ll + l2)))) % 32; + int gi3 = perm(ii + i3 + perm(jj + j3 + perm(kk + k3 + perm(ll + l3)))) % 32; + int gi4 = perm(ii + 1 + perm(jj + 1 + perm(kk + 1 + perm(ll + 1)))) % 32; + double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0; + if (t0 < 0) + n0 = 0.0; + else { + t0 *= t0; + n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0); + } + double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1; + if (t1 < 0) + n1 = 0.0; + else { + t1 *= t1; + n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1); + } + double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2; + if (t2 < 0) + n2 = 0.0; + else { + t2 *= t2; + n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2); + } + + double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3; + if (t3 < 0) + n3 = 0.0; + else { + t3 *= t3; + n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3); + } + double t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4; + if (t4 < 0) + n4 = 0.0; + else { + t4 *= t4; + n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4); + } + return 27.0 * (n0 + n1 + n2 + n3 + n4); + } + }; +} // namespace shammath diff --git a/src/shammath/src/Perlin.cpp b/src/shammath/src/Perlin.cpp new file mode 100755 index 0000000000..bb336dabc2 --- /dev/null +++ b/src/shammath/src/Perlin.cpp @@ -0,0 +1,288 @@ +/* + * perlin.cpp + * + * Created on: May 18, 2017 + * Author: tim + */ + + + + +#include "Perlin.h" + +namespace noise { + + int grad3[][3] = {{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0}, + {1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1}, + {0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}}; + + int grad4[][4]= {{0,1,1,1}, {0,1,1,-1}, {0,1,-1,1}, {0,1,-1,-1}, + {0,-1,1,1}, {0,-1,1,-1}, {0,-1,-1,1}, {0,-1,-1,-1}, + {1,0,1,1}, {1,0,1,-1}, {1,0,-1,1}, {1,0,-1,-1}, + {-1,0,1,1}, {-1,0,1,-1}, {-1,0,-1,1}, {-1,0,-1,-1}, + {1,1,0,1}, {1,1,0,-1}, {1,-1,0,1}, {1,-1,0,-1}, + {-1,1,0,1}, {-1,1,0,-1}, {-1,-1,0,1}, {-1,-1,0,-1}, + {1,1,1,0}, {1,1,-1,0}, {1,-1,1,0}, {1,-1,-1,0}, + {-1,1,1,0}, {-1,1,-1,0}, {-1,-1,1,0}, {-1,-1,-1,0}}; + + int p[] = {151,160,137,91,90,15, + 131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23, + 190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33, + 88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166, + 77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244, + 102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196, + 135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123, + 5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42, + 223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9, + 129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228, + 251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107, + 49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254, + 138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180}; + + int* perm = new int[512]; + + void init(){ + for(int i=0; i<512; i++){ + perm[i]=p[i & 255]; + } + } + + int simplex[][4] = { + {0,1,2,3},{0,1,3,2},{0,0,0,0},{0,2,3,1},{0,0,0,0},{0,0,0,0},{0,0,0,0},{1,2,3,0}, + {0,2,1,3},{0,0,0,0},{0,3,1,2},{0,3,2,1},{0,0,0,0},{0,0,0,0},{0,0,0,0},{1,3,2,0}, + {0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}, + {1,2,0,3},{0,0,0,0},{1,3,0,2},{0,0,0,0},{0,0,0,0},{0,0,0,0},{2,3,0,1},{2,3,1,0}, + {1,0,2,3},{1,0,3,2},{0,0,0,0},{0,0,0,0},{0,0,0,0},{2,0,3,1},{0,0,0,0},{2,1,3,0}, + {0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}, + {2,0,1,3},{0,0,0,0},{0,0,0,0},{0,0,0,0},{3,0,1,2},{3,0,2,1},{0,0,0,0},{3,1,2,0}, + {2,1,0,3},{0,0,0,0},{0,0,0,0},{0,0,0,0},{3,1,0,2},{0,0,0,0},{3,2,0,1},{3,2,1,0}}; + + + + + double noise(double xin, double yin){ + double n0, n1, n2; + + double F2 = 0.5*(sqrt(3.0)-1.0); + double s = (xin+yin)*F2; + + int i = fastfloor(xin+s); + int j = fastfloor(yin+s); + double G2 = (3.0-sqrt(3.0))/6.0; + double t = (i+j)*G2; + double X0 = i-t; + + double Y0 = j-t; + double x0 = xin-X0; + + double y0 = yin-Y0; + int i1, j1; + if(x0>y0) {i1=1; j1=0;} + else {i1=0; j1=1;} + double x1 = x0 - i1 + G2; + + double y1 = y0 - j1 + G2; + double x2 = x0 - 1.0 + 2.0 * G2; + double y2 = y0 - 1.0 + 2.0 * G2; + int ii = i & 255; + int jj = j & 255; + int gi0 = perm[ii+perm[jj]] % 12; + int gi1 = perm[ii+i1+perm[jj+j1]] % 12; + int gi2 = perm[ii+1+perm[jj+1]] % 12; + double t0 = 0.5 - x0*x0-y0*y0; + if(t0<0) n0 = 0.0; + else { + t0 *= t0; + n0 = t0 * t0 * dot(grad3[gi0], x0, y0); + } + double t1 = 0.5 - x1*x1-y1*y1; + if(t1<0) n1 = 0.0; + else { + t1 *= t1; + n1 = t1 * t1 * dot(grad3[gi1], x1, y1); + } + + double t2 = 0.5 - x2*x2-y2*y2; + if(t2<0) n2 = 0.0; + else { + t2 *= t2; + n2 = t2 * t2 * dot(grad3[gi2], x2, y2); + } + + return 70.0 * (n0 + n1 + n2); + } + + + double noise(double xin, double yin, double zin){ + double n0, n1, n2, n3; + double F3 = 1.0/3.0; + double s = (xin+yin+zin)*F3; + int i = fastfloor(xin+s); + int j = fastfloor(yin+s); + int k = fastfloor(zin+s); + double G3 = 1.0/6.0; + double t = (i+j+k)*G3; + double X0 = i-t; + double Y0 = j-t; + double Z0 = k-t; + double x0 = xin-X0; + double y0 = yin-Y0; + double z0 = zin-Z0; + int i1, j1, k1; + int i2, j2, k2; + if(x0>=y0) { + if(y0>=z0) + { i1=1; j1=0; k1=0; i2=1; j2=1; k2=0; } + else if(x0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=0; k2=1; } + else { i1=0; j1=0; k1=1; i2=1; j2=0; k2=1; } + } + else { + if(y0 y0) ? 32 : 0; + int c2 = (x0 > z0) ? 16 : 0; + int c3 = (y0 > z0) ? 8 : 0; + int c4 = (x0 > w0) ? 4 : 0; + int c5 = (y0 > w0) ? 2 : 0; + int c6 = (z0 > w0) ? 1 : 0; + int c = c1 + c2 + c3 + c4 + c5 + c6; + int i1, j1, k1, l1; + int i2, j2, k2, l2; + int i3, j3, k3, l3; + i1 = simplex[c][0]>=3 ? 1 : 0; + j1 = simplex[c][1]>=3 ? 1 : 0; + k1 = simplex[c][2]>=3 ? 1 : 0; + l1 = simplex[c][3]>=3 ? 1 : 0; + i2 = simplex[c][0]>=2 ? 1 : 0; + j2 = simplex[c][1]>=2 ? 1 : 0; + + k2 = simplex[c][2]>=2 ? 1 : 0; + l2 = simplex[c][3]>=2 ? 1 : 0; + i3 = simplex[c][0]>=1 ? 1 : 0; + j3 = simplex[c][1]>=1 ? 1 : 0; + k3 = simplex[c][2]>=1 ? 1 : 0; + l3 = simplex[c][3]>=1 ? 1 : 0; + double x1 = x0 - i1 + G4; + double y1 = y0 - j1 + G4; + double z1 = z0 - k1 + G4; + double w1 = w0 - l1 + G4; + double x2 = x0 - i2 + 2.0*G4; + double y2 = y0 - j2 + 2.0*G4; + double z2 = z0 - k2 + 2.0*G4; + double w2 = w0 - l2 + 2.0*G4; + double x3 = x0 - i3 + 3.0*G4; + double y3 = y0 - j3 + 3.0*G4; + double z3 = z0 - k3 + 3.0*G4; + double w3 = w0 - l3 + 3.0*G4; + double x4 = x0 - 1.0 + 4.0*G4; + double y4 = y0 - 1.0 + 4.0*G4; + double z4 = z0 - 1.0 + 4.0*G4; + double w4 = w0 - 1.0 + 4.0*G4; + int ii = i & 255; + int jj = j & 255; + int kk = k & 255; + int ll = l & 255; + int gi0 = perm[ii+perm[jj+perm[kk+perm[ll]]]] % 32; + int gi1 = perm[ii+i1+perm[jj+j1+perm[kk+k1+perm[ll+l1]]]] % 32; + int gi2 = perm[ii+i2+perm[jj+j2+perm[kk+k2+perm[ll+l2]]]] % 32; + int gi3 = perm[ii+i3+perm[jj+j3+perm[kk+k3+perm[ll+l3]]]] % 32; + int gi4 = perm[ii+1+perm[jj+1+perm[kk+1+perm[ll+1]]]] % 32; + double t0 = 0.6 - x0*x0 - y0*y0 - z0*z0 - w0*w0; + if(t0<0) n0 = 0.0; + else { + t0 *= t0; + n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0); + } + double t1 = 0.6 - x1*x1 - y1*y1 - z1*z1 - w1*w1; + if(t1<0) n1 = 0.0; + else { + t1 *= t1; + n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1); + } + double t2 = 0.6 - x2*x2 - y2*y2 - z2*z2 - w2*w2; + if(t2<0) n2 = 0.0; + else { + t2 *= t2; + n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2); + } + + double t3 = 0.6 - x3*x3 - y3*y3 - z3*z3 - w3*w3; + if(t3<0) n3 = 0.0; + else { + t3 *= t3; + n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3); + } + double t4 = 0.6 - x4*x4 - y4*y4 - z4*z4 - w4*w4; + if(t4<0) n4 = 0.0; + else { + t4 *= t4; + n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4); + } + return 27.0 * (n0 + n1 + n2 + n3 + n4); + } + +} diff --git a/src/shammodels/common/include/shammodels/common/EOSConfig.hpp b/src/shammodels/common/include/shammodels/common/EOSConfig.hpp index 1613b1fae1..c51bea0c57 100644 --- a/src/shammodels/common/include/shammodels/common/EOSConfig.hpp +++ b/src/shammodels/common/include/shammodels/common/EOSConfig.hpp @@ -66,6 +66,9 @@ namespace shammodels { using LocallyIsothermalFA2014 = shamphys::EOS_Config_LocallyIsothermalDisc_Farris2014; + /// Machida 2006 equation of state configuration + using Machida06 = shamphys::EOS_Config_Machida06; + /// Variant type to store the EOS configuration using Variant = std::variant< Isothermal, @@ -73,7 +76,8 @@ namespace shammodels { Polytropic, LocallyIsothermal, LocallyIsothermalLP07, - LocallyIsothermalFA2014>; + LocallyIsothermalFA2014, + Machida06>; /// Current EOS configuration Variant config = Adiabatic{}; @@ -123,6 +127,22 @@ namespace shammodels { config = LocallyIsothermalFA2014{h_over_r}; } + /** + * @brief Set the EOS configuration to a Machida 2006 equation of state + * + * @param rho_c1 Critical density 1 + * @param rho_c2 Critical density 2 + * @param rho_c3 Critical density 3 + * @param cs Sound speed + * @param mu Mean molecular weight + * @param mh Mass of hydrogen + * @param kb Boltzmann constant + */ + inline void set_machida06( + Tscal rho_c1, Tscal rho_c2, Tscal rho_c3, Tscal cs, Tscal mu, Tscal mh, Tscal kb) { + config = Machida06{rho_c1, rho_c2, rho_c3, cs, mu, mh, kb}; + } + /** * @brief Print current status of the EOSConfig */ diff --git a/src/shammodels/common/src/EOSConfig.cpp b/src/shammodels/common/src/EOSConfig.cpp index 962485c7fe..f16cccbeb2 100644 --- a/src/shammodels/common/src/EOSConfig.cpp +++ b/src/shammodels/common/src/EOSConfig.cpp @@ -61,7 +61,7 @@ namespace shammodels { using LocIsoT = typename EOSConfig::LocallyIsothermal; using LocIsoTLP07 = typename EOSConfig::LocallyIsothermalLP07; using LocIsoTFA2014 = typename EOSConfig::LocallyIsothermalFA2014; - + using Machida06 = typename EOSConfig::Machida06; if (const Isothermal *eos_config = std::get_if(&p.config)) { j = json{{"Tvec", type_id}, {"eos_type", "isothermal"}, {"cs", eos_config->cs}}; } else if (const Adiabatic *eos_config = std::get_if(&p.config)) { @@ -86,6 +86,17 @@ namespace shammodels { {"Tvec", type_id}, {"eos_type", "locally_isothermal_fa2014"}, {"h_over_r", eos_config->h_over_r}}; + } else if (const Machida06 *eos_config = std::get_if(&p.config)) { + j = json{ + {"Tvec", type_id}, + {"eos_type", "machida06"}, + {"rho_c1", eos_config->rho_c1}, + {"rho_c2", eos_config->rho_c2}, + {"rho_c3", eos_config->rho_c3}, + {"cs", eos_config->cs}, + {"mu", eos_config->mu}, + {"mh", eos_config->mh}, + {"kb", eos_config->kb}}; } else { shambase::throw_unimplemented(); // should never be reached } @@ -138,7 +149,7 @@ namespace shammodels { using LocIsoT = typename EOSConfig::LocallyIsothermal; using LocIsoTLP07 = typename EOSConfig::LocallyIsothermalLP07; using LocIsoTFA2014 = typename EOSConfig::LocallyIsothermalFA2014; - + using Machida06 = typename EOSConfig::Machida06; if (eos_type == "isothermal") { p.config = Isothermal{j.at("cs").get()}; } else if (eos_type == "adiabatic") { @@ -152,6 +163,15 @@ namespace shammodels { j.at("cs0").get(), j.at("q").get(), j.at("r0").get()}; } else if (eos_type == "locally_isothermal_fa2014") { p.config = LocIsoTFA2014{j.at("h_over_r").get()}; + } else if (eos_type == "machida06") { + p.config = Machida06{ + j.at("rho_c1").get(), + j.at("rho_c2").get(), + j.at("rho_c3").get(), + j.at("cs").get(), + j.at("mu").get(), + j.at("mh").get(), + j.at("kb").get()}; } else { shambase::throw_unimplemented("wtf !"); } diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp index fb9e98090d..7087b23284 100644 --- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp @@ -558,6 +558,22 @@ struct shammodels::sph::SolverConfig { eos_config.set_locally_isothermalFA2014(h_over_r); } + /** + * @brief Set the EOS configuration to a Machida 2006 equation of state + * + * @param rho_c1 Critical density 1 + * @param rho_c2 Critical density 2 + * @param rho_c3 Critical density 3 + * @param cs Sound speed + * @param mu Mean molecular weight + * @param mh Mass of hydrogen + * @param kb Boltzmann constant + */ + inline void set_eos_machida06( + Tscal rho_c1, Tscal rho_c2, Tscal rho_c3, Tscal cs, Tscal mu, Tscal mh, Tscal kb) { + eos_config.set_machida06(rho_c1, rho_c2, rho_c3, cs, mu, mh, kb); + } + ////////////////////////////////////////////////////////////////////////////////////////////// // EOS Config (END) ////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/shammodels/sph/src/modules/ComputeEos.cpp b/src/shammodels/sph/src/modules/ComputeEos.cpp index 8be298bd7a..4165e7903b 100644 --- a/src/shammodels/sph/src/modules/ComputeEos.cpp +++ b/src/shammodels/sph/src/modules/ComputeEos.cpp @@ -108,6 +108,7 @@ void shammodels::sph::modules::ComputeEos::compute_eos_internal using SolverEOS_LocallyIsothermal = typename SolverConfigEOS::LocallyIsothermal; using SolverEOS_LocallyIsothermalLP07 = typename SolverConfigEOS::LocallyIsothermalLP07; using SolverEOS_LocallyIsothermalFA2014 = typename SolverConfigEOS::LocallyIsothermalFA2014; + using SolverEOS_Machida06 = typename SolverConfigEOS::Machida06; sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); @@ -389,6 +390,58 @@ void shammodels::sph::modules::ComputeEos::compute_eos_internal buf_xyz.complete_event_state(e); }); + } else if ( + SolverEOS_Machida06 *eos_config + = std::get_if(&solver_config.eos_config.config)) { + + using EOS = shamphys::EOS_Machida06; + + storage.merged_patchdata_ghost.get().for_each([&](u64 id, PatchDataLayer &mpdat) { + auto &mfield = storage.merged_xyzh.get().get(id); + + sham::DeviceBuffer &buf_xyz = mfield.template get_field_buf_ref(0); + + sham::DeviceBuffer &buf_P + = shambase::get_check_ref(storage.pressure).get_field(id).get_buf(); + sham::DeviceBuffer &buf_cs + = shambase::get_check_ref(storage.soundspeed).get_field(id).get_buf(); + sham::DeviceBuffer &buf_uint = mpdat.get_field_buf_ref(iuint_interf); + auto rho_getter = rho_getter_gen(mpdat); + + Tscal rho_c1 = eos_config->rho_c1; + Tscal rho_c2 = eos_config->rho_c2; + Tscal rho_c3 = eos_config->rho_c3; + Tscal cs = eos_config->cs; + Tscal mu = eos_config->mu; + Tscal mh = eos_config->mh; + Tscal kb = eos_config->kb; + + u32 total_elements + = shambase::get_check_ref(storage.part_counts_with_ghost).indexes.get(id); + + sham::kernel_call( + q, + sham::MultiRef{rho_getter, buf_uint, buf_xyz}, + sham::MultiRef{buf_P, buf_cs}, + total_elements, + [rho_c1, rho_c2, rho_c3, cs0 = cs, mu, mh, kb]( + u32 i, + auto rho, + const Tscal *__restrict U, + const Tvec *__restrict xyz, + Tscal *__restrict P, + Tscal *__restrict cs) { + using namespace shamrock::sph; + + Tscal rho_a = rho(i); + + Tscal P_a = EOS::pressure(cs0, rho_a, rho_c1, rho_c2, rho_c3); + Tscal cs_a = EOS::soundspeed(P_a, rho_a, rho_c1, rho_c2, rho_c3); + + P[i] = P_a; + cs[i] = cs_a; + }); + }); } else { shambase::throw_unimplemented(); } diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index b63441bc8c..38311ddca9 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -95,6 +95,26 @@ void add_instance(py::module &m, std::string name_config, std::string name_model }, py::kw_only(), py::arg("h_over_r")) + .def( + "set_eos_machida06", + [](TConfig &self, + Tscal rho_c1, + Tscal rho_c2, + Tscal rho_c3, + Tscal cs, + Tscal mu, + Tscal mh, + Tscal kb) { + self.set_eos_machida06(rho_c1, rho_c2, rho_c3, cs, mu, mh, kb); + }, + py::kw_only(), + py::arg("rho_c1"), + py::arg("rho_c2"), + py::arg("rho_c3"), + py::arg("cs"), + py::arg("mu"), + py::arg("mh"), + py::arg("kb")) .def("set_artif_viscosity_None", &TConfig::set_artif_viscosity_None) .def( "to_json", diff --git a/src/shamphys/include/shamphys/eos_config.hpp b/src/shamphys/include/shamphys/eos_config.hpp index 7016a73632..6d15772ec4 100644 --- a/src/shamphys/include/shamphys/eos_config.hpp +++ b/src/shamphys/include/shamphys/eos_config.hpp @@ -187,4 +187,41 @@ namespace shamphys { return (lhs.h_over_r == rhs.h_over_r); } + /** + * @brief Configuration struct for the Machida 2006 equation of state + * + * @tparam Tscal Scalar type + */ + template + struct EOS_Config_Machida06 { + Tscal rho_c1; + Tscal rho_c2; + Tscal rho_c3; + Tscal cs; + Tscal mu; + Tscal mh; + Tscal kb; + }; + + /** + * @brief Equal operator for the EOS_Config_Machida06 struct + * + * @tparam Tscal Scalar type + * @param lhs First EOS_Config_Machida06 struct to compare + * @param rhs Second EOS_Config_Machida06 struct to compare + * + * This function checks if two EOS_Config_Machida06 structs are equal by comparing their rho_c1, + * rho_c2, rho_c3, cs, mu, mh, and kb values. + * + * @return true if the two structs have the same rho_c1, rho_c2, rho_c3, cs, mu, mh, and kb + * values, false otherwise + */ + template + inline bool operator==( + const EOS_Config_Machida06 &lhs, const EOS_Config_Machida06 &rhs) { + return (lhs.rho_c1 == rhs.rho_c1) && (lhs.rho_c2 == rhs.rho_c2) + && (lhs.rho_c3 == rhs.rho_c3) && (lhs.cs == rhs.cs) && (lhs.mu == rhs.mu) + && (lhs.mh == rhs.mh) && (lhs.kb == rhs.kb); + } + } // namespace shamphys diff --git a/src/shampylib/src/math/pyshammath.cpp b/src/shampylib/src/math/pyshammath.cpp index f108531bbf..760457809f 100644 --- a/src/shampylib/src/math/pyshammath.cpp +++ b/src/shampylib/src/math/pyshammath.cpp @@ -18,6 +18,7 @@ #include "shambackends/typeAliasVec.hpp" #include "shambindings/pybindaliases.hpp" #include "shambindings/pytypealias.hpp" +#include "shammath/Perlin.h" #include "shammath/derivatives.hpp" #include "shammath/matrix.hpp" #include "shammath/matrix_op.hpp" @@ -849,4 +850,9 @@ Register_pymod(pysham_mathinit) { "SymTensorCollection_f64_1_1(\n t1={}\n)", py::str(py::cast(c.t1)).cast()); }); + + // PerlinNoise + py::class_(math_module, "PerlinNoise") + .def(py::init<>()) + .def("noise_3d", &shammath::PerlinNoise::noise_3d); }