yohai/phigrape
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Improved printing of external potential

Browse source
This commit is contained in:
Yohai Meiron 2020-10-31 15:08:48 -04:00
parent 8adb4ac813
commit dd31b63215
3 changed files with 57 additions and 35 deletions
Download patch file Download diff file Expand all files Collapse all files

2
TODO.md
View file

@ -4,3 +4,5 @@ TODO
* Break main() into smaller chunks; operations that are timed should become independent functions.
* Const everything
* OpenMP

39
external.h
View file

@ -1,4 +1,5 @@
#pragma once
#include <string>
#include <vector>
#include "double3.h"
@ -15,11 +16,17 @@ public:
}
}
virtual void calc_gravity() = 0;
bool is_active;
virtual void print_info() {}
void set_name(const std::string &name)
{
this->name = name;
}
bool is_active = false;
protected:
double potential;
double3 acceleration, jerk;
double3 x, v;
std::string name = "ext";
void set_coordinates(double3 x, double3 v)
{
this->x = x;
@ -30,23 +37,38 @@ protected:
class Plummer : public External_gravity {
public:
Plummer(double m, double b) : m(m), b(b) {is_active=(m>0);}
void calc_gravity();
void calc_gravity() override;
void print_info() override
{
if (!is_active) return;
printf("m_%-5s = %.4E b_%-5s = %.4E\n", name.c_str(), m, name.c_str(), b);
}
private:
double m, b;
};
class Miyamoto_Nagai : public External_gravity {
public:
Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0);}
Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0); this->set_name("disk");}
void calc_gravity();
void print_info() override
{
if (!is_active) return;
printf("m_%-5s = %.4E a_%-5s = %.4E b_%-5s = %.4E\n", name.c_str(), m, name.c_str(), a, name.c_str(), b);
}
private:
double m, a, b;
};
class Logarithmic_halo : public External_gravity {
public:
Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0);}
void calc_gravity();
Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0); this->set_name("halo");}
void calc_gravity() override;
void print_info() override
{
if (!is_active) return;
printf("v_%-4s = %.6E r_%-4s = %.4E\n", name.c_str(), sqrt(v2_halo), name.c_str(), sqrt(r2_halo));
}
private:
double v2_halo, r2_halo;
};
@ -54,7 +76,12 @@ private:
class Dehnen : public External_gravity {
public:
Dehnen(double m, double r, double gamma) : m(m), r(r), gamma(gamma) {is_active=(m>0);}
void calc_gravity();
void calc_gravity() override;
void print_info() override
{
if (!is_active) return;
printf("m_%-5s = %.4E r_%-5s = %.4E g_%-5s = %.4E\n", name.c_str(), m, name.c_str(), r, name.c_str(), gamma);
}
private:
double m, r, gamma;
};

43
phigrape.cpp
View file

@ -85,6 +85,13 @@ public:
component->apply(n, x, v, pot, acc, jrk);
}
}
void print_info()
{
for (auto component : components) {
component->print_info();
}
fflush(stdout);
}
private:
std::vector<External_gravity*> components;
};
@ -249,13 +256,8 @@ inline double aarseth_step(const double eta, const double dt, const double3 &a,
int main(int argc, char *argv[])
{
/* read the input parameters */
const Config config("phigrape.conf");
timer.start();
double timesteps=0.0, n_act_sum=0.0;
/* INIT the rand() !!! */
srand(19640916); /* it is just my birthday :-) */
@ -276,6 +278,7 @@ int main(int argc, char *argv[])
/* Print the Rank and the names of processors */
printf("Rank of the processor %03d on %s \n", myRank, processor_name);
const Config config("phigrape.conf");
Input_data input_data;
if (is_hdf5(config.input_file_name)) {
#ifndef HAS_HDF5
@ -294,6 +297,10 @@ int main(int argc, char *argv[])
auto &x = input_data.x;
auto &v = input_data.v;
double t_disk = config.dt_disk*(1+floor(time_cur/config.dt_disk));
double t_contr = config.dt_contr*(1+floor(time_cur/config.dt_contr));
double t_bh = config.dt_bh*(1+floor(time_cur/config.dt_bh));
if (myRank == rootRank) {
printf("\n");
printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc);
@ -302,8 +309,8 @@ int main(int argc, char *argv[])
printf("t_beg = %.6E \t t_end = %.6E\n", time_cur, config.t_end);
printf("dt_disk = %.6E \t dt_contr = %.6E\n", config.dt_disk, config.dt_contr);
printf("dt_bh = %.6E \n", config.dt_bh);
printf("eta = %.6E \n", config.eta);
printf("\n");
printf("eta = %.6E\n\n", config.eta);
printf("t_disk = %.6E t_contr = %.6E t_bh = %.6E\n\n", t_disk, t_contr, t_bh);
if ((diskstep == 0) && (time_cur == 0)) {
FILE *out = fopen("contr.dat", "w");
@ -327,33 +334,18 @@ int main(int argc, char *argv[])
}
Calc_ext_grav calc_ext_grav;
Plummer ext_bulge(config.ext_m_bulge*normalization_mass, config.ext_b_bulge*normalization_length);
ext_bulge.set_name("bulge");
calc_ext_grav.add_component(ext_bulge);
Miyamoto_Nagai ext_disk(config.ext_m_disk*normalization_mass, config.ext_a_disk*normalization_length, config.ext_b_disk*normalization_length);
calc_ext_grav.add_component(ext_disk);
Plummer ext_halo_plummer(config.ext_m_halo_plummer*normalization_mass, config.ext_b_halo_plummer*normalization_length);
ext_halo_plummer.set_name("halo");
calc_ext_grav.add_component(ext_halo_plummer);
Logarithmic_halo ext_log_halo(config.ext_log_halo_v*normalization_velocity, config.ext_log_halo_r*normalization_length);
calc_ext_grav.add_component(ext_log_halo);
Dehnen ext_dehnen(config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma);
calc_ext_grav.add_component(ext_dehnen);
bool has_external_gravity = calc_ext_grav.any_active();
if (has_external_gravity) printf("External Potential: \n\n");
if (ext_bulge.is_active) printf("m_bulge = %.4E b_bulge = %.4E\n", config.ext_m_bulge*normalization_mass, config.ext_b_bulge*normalization_length);
if (ext_disk.is_active) printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", config.ext_m_disk*normalization_mass, config.ext_a_disk*normalization_length, config.ext_b_disk*normalization_length);
if (ext_halo_plummer.is_active) printf("m_halo = %.4E b_halo = %.4E\n", config.ext_m_halo_plummer*normalization_mass, config.ext_b_halo_plummer*normalization_length);
if (ext_log_halo.is_active) printf("v_halo = %.6E r_halo = %.6E \n", config.ext_log_halo_v*normalization_velocity, config.ext_log_halo_r*normalization_length);
if (ext_dehnen.is_active) printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n\n", config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma);
fflush(stdout);
double t_disk = config.dt_disk*(1+floor(time_cur/config.dt_disk));
double t_contr = config.dt_contr*(1+floor(time_cur/config.dt_contr));
double t_bh = config.dt_bh*(1+floor(time_cur/config.dt_bh));
if (myRank == rootRank) {
printf("t_disk = %.6E t_contr = %.6E t_bh = %.6E \n\n", t_disk, t_contr, t_bh);
fflush(stdout);
} /* if (myRank == rootRank) */
if (myRank == rootRank) calc_ext_grav.print_info();
/* some local settings for G6a boards */
int clusterid, numGPU;
@ -439,6 +431,7 @@ int main(int argc, char *argv[])
std::vector<double> pot_ext(N, 0.);
calc_ext_grav(N, x, v, pot_ext, a, adot);
double timesteps=0, n_act_sum=0;
/* Energy control... */
if (myRank == rootRank) energy_contr(time_cur, timesteps, n_act_sum, calc_self_grav.g6_calls, N, m, x, v, pot, pot_ext);