yohai/etics-external
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Initial commit

Browse source
This commit is contained in:
Yohai Meiron 2025-11-04 22:25:38 -05:00
commit b59b6716a6
3 changed files with 173 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

5
Makefile Normal file
View file

@ -0,0 +1,5 @@
main:
$(CXX) -O3 main.cpp ../grapite/libgrapite.a -lcuda -lcudart -o $@
clean:
$(RM) main

58
finish_processing.py Normal file
View file

@ -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()

110
main.cpp Normal file
View file

@ -0,0 +1,110 @@
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
#include "../grapite/etics_interface.h"
#include "../grapite/particle_manager.h"
// #include <format> // C++23
//#include <print> // C++23
#include <cstdio> // 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<long long> 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<double>(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<double> pot(n_total);
std::vector<vec3> 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<double>(counter)*100);
fflush(stdout);
}
}
printf("\rFinished\n");
}