Skip to content

Commit

Permalink
Merge pull request #1086 from LourensVeen/issue-940-update-mpreal-h
Browse files Browse the repository at this point in the history 
Issue 940: update mpreal.h
  • Loading branch information
@LourensVeen
LourensVeen authored Oct 31, 2024
2 parents 8f44615 + 4b6b016 commit 411ca13
Show file tree
Hide file tree
Showing 12 changed files with 7,517 additions and 6,477 deletions.
6,570 changes: 3,377 additions & 3,193 deletions src/amuse/community/adaptb/src/mpreal.h
View file

Large diffs are not rendered by default.

62 changes: 32 additions & 30 deletions src/amuse/community/brutus/interface.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,6 @@
// Includes
////////////////////////////////////////////////////////
#include <iostream>
using namespace std;

#include <fstream>
#include <cstdlib>
Expand All @@ -14,21 +13,21 @@ using namespace std;
////////////////////////////////////////////////////////
// Declare global variables
////////////////////////////////////////////////////////
ofstream odata;
std::ofstream odata;

int particle_id_counter = 0;

int numBits = 64;
int numBits = 64;
int numDigits = numBits/4;

string out_directory;
std::string out_directory;
std::map<int, int> local_index_map;

/*
* We need this result_strings array to ensure that
* C++ strings are not reclaimed before the function ends
*/
string result_strings[10];
std::string result_strings[10];

Brutus *brutus = NULL;

Expand All @@ -38,15 +37,16 @@ mpreal t = "0";

mpreal epsilon = "1e-6"; // Bulirsch-Stoer tolerance

vector<mpreal> data, data_radius;
std::vector<mpreal> data, data_radius;

////////////////////////////////////////////////////////
// Amuse interface functions
////////////////////////////////////////////////////////

int initialize_code() {
odata.open("temp.log");

mpreal::set_default_prec(numBits);
mpreal::set_default_prec(numBits);

brutus = new Brutus();

Expand All @@ -62,7 +62,7 @@ int initialize_code() {

// functions with "_string" assign strings to mpreals, and without "_string" assign doubles to mpreals

int new_particle_string(int *particle_identifier, char* mass,
int new_particle_string(int *particle_identifier, char* mass,
char* x, char* y, char* z, char* vx, char* vy, char* vz, char* radius) {

data.push_back(mass);
Expand All @@ -80,7 +80,7 @@ int new_particle_string(int *particle_identifier, char* mass,

return 0;
}
int new_particle_float64(int *particle_identifier, double mass,
int new_particle_float64(int *particle_identifier, double mass,
double x, double y, double z, double vx, double vy, double vz, double radius) {

data.push_back( (mpreal)mass );
Expand Down Expand Up @@ -139,7 +139,7 @@ int get_begin_time(double * output) {
// Timestep parameter, eta
int set_eta_string(char* myeta) {
eta = myeta;
brutus->set_eta(eta);
brutus->set_eta(eta);
return 0;
}
int get_eta_string(char **myeta) {
Expand Down Expand Up @@ -200,12 +200,12 @@ int get_bs_tolerance_string(char **bs_tolerance) {
}
int set_bs_tolerance(double bs_tolerance) {

odata << t << ": changing e from " << epsilon << ", to " << bs_tolerance << endl;
odata << t << ": changing e from " << epsilon << ", to " << bs_tolerance << std::endl;

epsilon = (mpreal) bs_tolerance;
brutus->set_tolerance(epsilon);

odata << epsilon << " " << brutus->get_tolerance() << endl;
odata << epsilon << " " << brutus->get_tolerance() << std::endl;

return 0;
}
Expand All @@ -215,14 +215,14 @@ int get_bs_tolerance(double *bs_tolerance) {
}
// Word-length, numBits in mantissa
int set_word_length(int mynumBits) {
odata << t << ": changing L from " << numBits << ", to " << mynumBits << endl;
odata << t << ": changing L from " << numBits << ", to " << mynumBits << std::endl;

numBits = mynumBits;
mpreal::set_default_prec(numBits);
numDigits = (int)abs(log10( pow("2.0", -numBits) )).toLong();
brutus->set_numBits(numBits);

odata << numBits << " " << brutus->get_numBits() << endl;
odata << numBits << " " << brutus->get_numBits() << std::endl;

return 0;
}
Expand Down Expand Up @@ -263,7 +263,7 @@ int get_mass_string(int id, char **mass) {
mass_string = data[id*7+0].toString();
*mass = (char*) mass_string.c_str();
return 0;
}
}
int set_mass_string(int id, char *mass) {
if (id < 0 || id >= particle_id_counter){
return -1;
Expand All @@ -277,7 +277,7 @@ int get_mass(int id, double* mass) {
}
*mass = data[id*7+0].toDouble();
return 0;
}
}
int set_mass(int id, double mass) {
if (id < 0 || id >= particle_id_counter){
return -1;
Expand Down Expand Up @@ -379,7 +379,7 @@ std::string get_state_strings_vx;
std::string get_state_strings_vy;
std::string get_state_strings_vz;
std::string get_state_strings_r;
int get_state_string(int id, char** m, char** x, char** y, char** z, char** vx, char** vy, char** vz, char** radius) {
int get_state_string(int id, char** m, char** x, char** y, char** z, char** vx, char** vy, char** vz, char** radius) {
if (id < 0 || id >= particle_id_counter){
return -1;
}
Expand Down Expand Up @@ -445,7 +445,7 @@ int set_state(int id, double m, double x, double y, double z, double vx, double
}

std::string radius_string;
int get_radius_string(int id, char** radius){
int get_radius_string(int id, char** radius){
if (id < 0 || id >= particle_id_counter){
return -1;
}
Expand All @@ -460,7 +460,7 @@ int set_radius_string(int id, char* radius) {
data_radius[id] = radius;
return 0;
}
int get_radius(int id, double* radius){
int get_radius(int id, double* radius){
if (id < 0 || id >= particle_id_counter){
return -1;
}
Expand Down Expand Up @@ -523,7 +523,7 @@ int get_index_of_first_particle(int* id){return -2;}
int get_index_of_next_particle(int id, int* idnext){return -2;}

std::string total_mass_string;
int get_total_mass_string(char **M){
int get_total_mass_string(char **M){
int N = data.size()/7;
mpreal Mtot = "0";
for(int i=0; i<N; i++) {
Expand All @@ -534,7 +534,7 @@ int get_total_mass_string(char **M){
return 0;
}

int get_total_mass(double* M){
int get_total_mass(double* M){
int N = data.size()/7;
mpreal Mtot = "0";
for(int i=0; i<N; i++) {
Expand Down Expand Up @@ -564,7 +564,7 @@ int get_potential_energy_m(mpreal* ep) {

eptot -= mi*mj/sqrt(dr2);
}
}
}

*ep = eptot;
return 0;
Expand Down Expand Up @@ -608,8 +608,8 @@ int get_kinetic_energy_string( char **ep) {
std::string total_energy_string;

int get_total_energy_string( char **ep) {
mpreal ektot = "0";
mpreal eptot = "0";
mpreal ektot = "0";
mpreal eptot = "0";
mpreal etot = "0";
get_potential_energy_m(&eptot);
get_kinetic_energy_m(&ektot);
Expand All @@ -633,40 +633,40 @@ int get_kinetic_energy(double* ek) {
return 0;
}

int get_center_of_mass_position(double* x , double* y, double* z){
int get_center_of_mass_position(double* x , double* y, double* z){
int N = data.size()/7;
mpreal Mtot = "0";
for(int i=0; i<N; i++) {
Mtot += data[i*7];
}

vector<mpreal> rcm(3,"0");
std::vector<mpreal> rcm(3,"0");
for(int i=0; i<N; i++) {
for(int j=0; j<3; j++) {
rcm[j] += data[i*7]*data[i*7+(j+1)];
}
}
for(int i=0; i<3; i++) rcm[i] /= Mtot;
for(int i=0; i<3; i++) rcm[i] /= Mtot;

*x = rcm[0].toDouble();
*y = rcm[1].toDouble();
*z = rcm[2].toDouble();
return 0;
}
int get_center_of_mass_velocity(double* vx, double* vy, double* vz){
int get_center_of_mass_velocity(double* vx, double* vy, double* vz){
int N = data.size()/7;
mpreal Mtot = "0";
for(int i=0; i<N; i++) {
Mtot += data[i*7];
}

vector<mpreal> vcm(3,"0");
std::vector<mpreal> vcm(3,"0");
for(int i=0; i<N; i++) {
for(int j=0; j<3; j++) {
vcm[j] += data[i*7]*data[i*7+(j+4)];
}
}
for(int i=0; i<3; i++) vcm[i] /= Mtot;
for(int i=0; i<3; i++) vcm[i] /= Mtot;

*vx = vcm[0].toDouble();
*vy = vcm[1].toDouble();
Expand All @@ -675,3 +675,5 @@ int get_center_of_mass_velocity(double* vx, double* vy, double* vz){
}

int get_acceleration(int id, double* ax, double* ay, double* az){return -2;}
int set_acceleration(int id, double ax, double ay, double az){return -2;}

Loading

0 comments on commit 411ca13

Please sign in to comment.