commit b59b6716a66889cb368d3740b791108b569cf148 Author: Yohai Meiron Date: Tue Nov 4 22:25:38 2025 -0500 Initial commit diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..0886a1b --- /dev/null +++ b/Makefile @@ -0,0 +1,5 @@ +main: + $(CXX) -O3 main.cpp ../grapite/libgrapite.a -lcuda -lcudart -o $@ + +clean: + $(RM) main \ No newline at end of file diff --git a/finish_processing.py b/finish_processing.py new file mode 100644 index 0000000..1a23062 --- /dev/null +++ b/finish_processing.py @@ -0,0 +1,58 @@ +#!/usr/bin/env python +import sys +import numpy as np, pandas as pd + +import importlib.util +spec = importlib.util.spec_from_file_location('scf', '../etics/src/scf.py') +scf = importlib.util.module_from_spec(spec) +spec.loader.exec_module(scf) + +n_max = 15 +l_max = 2 +length_scale = 10 + +pos_filename = '../../../sobolenko/TNG50-copy/for_yohai/98990136-weight.dat' +grav_filename = 'out.dat' +coeff_filename = 'out.coeff' + +pos_data = np.genfromtxt(pos_filename, skip_header=1, names=['idx', 'm', 'x', 'y', 'z', 'vx', 'vy', 'vz', 'eps', 'type', 'name', 't'], dtype=[int, float, float, float, float, float, float, float, float, int, int, int]) +grav_data = np.genfromtxt(grav_filename, names=['name', 'pot', 'ax', 'ay', 'az'], dtype=[int, float, float, float, float]) + +coeff_count = (n_max+1)*(l_max+1)*(l_max+2)//2 +A = np.fromfile(coeff_filename, dtype=np.complex128)[:coeff_count] +expansion = scf.Scf(n_max, l_max, length_scale=length_scale) +expansion.set_coefficients(A) + +r = np.sqrt(pos_data['x']**2 + pos_data['y']**2 + pos_data['z']**2) +r[r==0] = 1 # Will be ignored later +density_all = expansion.calc_radial_density(r) + +n = len(pos_data['idx']) +i2 = 0 +out_file = open('results.dat', 'w') +for i1 in range(10): + template = '%08d%15.8e%16.8e%16.8e%16.8e%16.8e%16.8e%16.8e%11.4e %d %015d%5.04d%16.8e%16.8e%16.8e%16.8e%15.8e\n' + if pos_data['name'][i1] == grav_data['name'][i2]: + pot, ax, ay, az = grav_data['pot'][i2], grav_data['ax'][i2], grav_data['ay'][i2], grav_data['az'][i2] + den = density_all[i1] + i2 += 1 + else: # Missing particle + pot, ax, ay, az = 0, 0, 0, 0 + den = 0 + out = template % (pos_data['idx'][i1], + pos_data['m'][i1], + pos_data['x'][i1], + pos_data['y'][i1], + pos_data['z'][i1], + pos_data['vx'][i1], + pos_data['vy'][i1], + pos_data['vz'][i1], + pos_data['eps'][i1], + pos_data['type'][i1], + pos_data['name'][i1], + pos_data['t'][i1], + pot, ax, ay, az, den + ) + out_file.write(out) + +out_file.close() \ No newline at end of file diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..40a84fb --- /dev/null +++ b/main.cpp @@ -0,0 +1,110 @@ +#include +#include +#include +#include +#include "../grapite/etics_interface.h" +#include "../grapite/particle_manager.h" + +// #include // C++23 +//#include // C++23 +#include // Legacy printf + +int main(int argc, char* argv[]) +{ + double _double; + int _int; + + if (argc != 2) { + std::cout << "Specify the input file\n"; + exit(1); + } + std::fstream file{argv[1], std::ios::in}; + if (!file) { + throw std::runtime_error("Could not open file"); + } + + int n_total; + std::string line{}; + std::getline(file, line); + std::stringstream s{line}; + s >> n_total >> _int >> _int >> _int; + if (!s) ("Error in header"); + printf("n_total = %d\n", n_total); + std::vector names; + + grapite::Passiv_particles particles_j{n_total}; + grapite::Active_particles particles_i{n_total}; // This is just needed to caluclate the forces, not the coefficients + int counter{}; + printf("Loading particles...\n"); + while (std::getline(file, line)) { + if (counter >= n_total) throw std::runtime_error("Particle exceeded maximum"); + //std::stringstream s{line}; + s = std::stringstream{line}; + int id; + double m, x, y, z; + long long name; + std::string _string; + s >> id >> m >> x >> y >> z >> _double >> _double >> _double >> _double >> _int >> name >> _int; + if (!s) throw std::runtime_error("Error in input"); + if (id == 0) continue; // skip the black hole + if (x == 0) continue; // bad particle + grapite::Particle_ext p{}; + p.m = m; + p.t = 0; + p.pos = vec3(x, y, z); + + particles_j.add_particle(counter, &p); + particles_i.add_particle(counter, &p); + names.emplace_back(name); + counter++; + if ((n_total > 100) && (counter % (n_total/100) == 0)) { + printf("\r%.0f%%", counter/static_cast(n_total)*100); + fflush(stdout); + } + } + n_total = counter; + printf("\rFinished\n"); + + // Initialize ETICS + grapite::Etics_interface etics{n_total}; + printf("ETICS was compiled with n_max=%d and l_max=%d\n", etics.GetNMAX(), etics.GetLMAX()); + constexpr double length_scale = 10.0; + etics::scf::SetLengthScale(length_scale); + printf("ETICS scaling is %.3e\n", length_scale); + + // Calculate the expanstion from particles_j + particles_j.nj_total = n_total; + particles_j.nj_core = 0; + particles_j.flush_memory(); + etics.calculate_expansion(&particles_j, grapite::Selector::halo, 0); + + // Write coefficients to a file + std::ofstream coeff_file("out.coeff", std::ofstream::binary); + coeff_file.write((char*)etics.A_h, 1024*sizeof(Complex)); + coeff_file.close(); + + // Calculate potentials and forces + particles_i.ni_core = 0; + particles_i.ni_total = n_total; + particles_i.flush_memory(); + std::vector pot(n_total); + std::vector acc(n_total); + etics.calculate_gravity(&particles_i, grapite::Selector::halo, 0, pot.data(), acc.data()); + + // Write results to a file + printf("Writing results...\n"); + std::ofstream output_file("out.dat"); + for (int i = 0; i < counter; i++) { + constexpr int buf_size = 4096; + char line[buf_size]; + auto [ax, ay, az] = acc[i]; + snprintf(line, buf_size, "%015lld%16.8e%16.8e%16.8e%16.8e\n", names[i], pot[i], ax, ay, az); + output_file.write(line, strlen(line)); + if ((counter > 100) && (i % (counter/100) == 0)) { + printf("\r%.0f%%", i/static_cast(counter)*100); + fflush(stdout); + } + } + printf("\rFinished\n"); + +} \ No newline at end of file