diff --git a/src/Makefile b/src/Makefile index 96bde0d..0379eee 100644 --- a/src/Makefile +++ b/src/Makefile @@ -2,7 +2,7 @@ PS_PATH = -I ../FDPS/src/ CC = g++ #CC = g++-mp-6 #CC = mpicxx-openmpi-gcc6 -CFLAGS = -O3 -ffast-math -funroll-loops +CFLAGS = -O3 -ffast-math -funroll-loops -fpermissive CFLAGS += -DPARTICLE_SIMULATOR_THREAD_PARALLEL -fopenmp #CFLAGS += -DPARTICLE_SIMULATOR_MPI_PARALLEL #CFLAGS += -Wall -Wformat=2 -Wcast-qual -Wcast-align -Wwrite-strings -Wconversion -Wfloat-equal -Wpointer-arith diff --git a/src/class.h b/src/class.h index b0aaa86..80f3501 100644 --- a/src/class.h +++ b/src/class.h @@ -27,6 +27,8 @@ class FileHeader{ } }; +// This namespace name can be confusing, because it is very close to std:: +// What about FDPS_SPH? Or whatever name you decide on. namespace STD{ namespace RESULT{ //Density summation @@ -78,6 +80,13 @@ namespace STD{ class RealPtcl{ public: + + // Instead of using abbreviations and explaining them in the comments, + // we have found it is much easier to read the code if you just use the + // full names. 'dens' is not much shorter than 'density', but 'density' + // is much easier to read (the same holds true for all other variables below). + // If I would encounter 'snds' in the code I would need to look up what it means + // but 'sound_speed' would be self-explanatory. PS::F64 mass; PS::F64vec pos, vel, acc; PS::F64 dens;//DENSity @@ -389,6 +398,16 @@ namespace STD{ } } +// As far as I can see the only 'Problem' that is currently implemented is init/GI.h? +// This class acts as the base class for any problem that might be implemented later, +// and it technically acts as something that is called an 'Interface' class. See +// https://www.tutorialspoint.com/cplusplus/cpp_interfaces.htm for a description of +// interface classes. Thus, I would recommend +// moving it to a new file called init/interface.h, and make it the only content +// of this file. +// Additionally, I find 'init' might not be fully explaining the content of that folder, it is not +// only describing the initial state, but also the equation of state, and the forces (and the domain). +// What about renaming it to 'model_setup' or 'problem_setup'? template class Problem{ Problem(){ } diff --git a/src/init/GI.h b/src/init/GI.h index 65423a4..5d96032 100644 --- a/src/init/GI.h +++ b/src/init/GI.h @@ -1,4 +1,3 @@ -#define SELF_GRAVITY #define FLAG_GI #ifdef PARTICLE_SIMULATOR_TWO_DIMENSION #error @@ -6,6 +5,8 @@ template class GI : public Problem{ public: static const double END_TIME = 1.0e+4; + static const bool self_gravity = true; + static void setupIC(PS::ParticleSystem& sph_system, system_t& sysinfo, PS::DomainInfo& dinfo){ const bool createTarget = true;//set false if you make an impactor. const double Corr = .98;//Correction Term diff --git a/src/main.cpp b/src/main.cpp index 2a79048..0dffa91 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include "param.h" #include "mathfunc.h" @@ -18,6 +19,17 @@ int main(int argc, char* argv[]){ namespace PTCL = STD; + + // There should be a capitalization rule for names, you currently use ALLCAPS + // for namespaces, but here you also use it for GI which is a class type. I would + // recommend using 'CamelCase' for everything except variable names, which you already + // name as 'names_with_underscores', but that is a personal preference and you are welcome + // to use your own. + + // Also there needs to be a way to include tests into the program. For example you could add a + // problem description called 'test' instead of 'GI' and if the program is started with a '--test' + // flag this test setup is executed. The results could then be compared to a reference test result. + typedef GI PROBLEM; ////////////////// //Create vars. @@ -56,27 +68,37 @@ int main(int argc, char* argv[]){ PS::TreeForForceShort::Gather dens_tree; PS::TreeForForceShort::Gather drvt_tree; PS::TreeForForceShort::Symmetry hydr_tree; - #ifdef SELF_GRAVITY - PS::TreeForForceLong ::Monopole grav_tree; - #endif + + // It would be advantageous if you could settle on a single style for input parameters, currently + // the code mixes variables, templates, and preprocessor macros. Ultimately it would be nice to + // have text files outside the source code that specify the model, but settling on a single header + // file would also be fine for the moment. If you want to save memory by not allocating a tree if this + // parameter is not set, consider making it an owning pointer (std::unique_ptr) like the following: + + std::unique_ptr::Monopole> grav_tree; + if (PROBLEM::self_gravity) + grav_tree.reset(new PS::TreeForForceLong ::Monopole()); + + // This way only the memory for a single pointer is allocated if the tree is not needed, and you can + // check for that by using if (grav_tree) dens_tree.initialize(sph_system.getNumberOfParticleLocal(), 0.5, 4, 128); drvt_tree.initialize(sph_system.getNumberOfParticleLocal(), 0.5, 4, 128); hydr_tree.initialize(sph_system.getNumberOfParticleLocal(), 0.5, 4, 128); - #ifdef SELF_GRAVITY - grav_tree.initialize(sph_system.getNumberOfParticleLocal(), 0.5, 8, 256); - #endif - for(short int loop = 0 ; loop <= PARAM::NUMBER_OF_DENSITY_SMOOTHING_LENGTH_LOOP ; ++ loop){ + if (PROBLEM::self_gravity) + grav_tree->initialize(sph_system.getNumberOfParticleLocal(), 0.5, 8, 256); + + for(short int loop = 0 ; loop <= PARAM::NUMBER_OF_DENSITY_SMOOTHING_LENGTH_LOOP ; ++ loop){ dens_tree.calcForceAllAndWriteBack(PTCL::CalcDensity(), sph_system, dinfo); } PTCL::CalcPressure(sph_system); drvt_tree.calcForceAllAndWriteBack(PTCL::CalcDerivative(), sph_system, dinfo); hydr_tree.calcForceAllAndWriteBack(PTCL::CalcHydroForce(), sph_system, dinfo); - #ifdef SELF_GRAVITY - grav_tree.calcForceAllAndWriteBack(PTCL::CalcGravityForce(), PTCL::CalcGravityForce(), sph_system, dinfo); - #endif - sysinfo.dt = getTimeStepGlobal(sph_system); + if (PROBLEM::self_gravity) + grav_tree->calcForceAllAndWriteBack(PTCL::CalcGravityForce(), PTCL::CalcGravityForce(), sph_system, dinfo); + + sysinfo.dt = getTimeStepGlobal(sph_system); PROBLEM::addExternalForce(sph_system, sysinfo); OutputFileWithTimeInterval(sph_system, sysinfo, PROBLEM::END_TIME); @@ -106,10 +128,10 @@ int main(int argc, char* argv[]){ PTCL::CalcPressure(sph_system); drvt_tree.calcForceAllAndWriteBack(PTCL::CalcDerivative(), sph_system, dinfo); hydr_tree.calcForceAllAndWriteBack(PTCL::CalcHydroForce(), sph_system, dinfo); - #ifdef SELF_GRAVITY - grav_tree.calcForceAllAndWriteBack(PTCL::CalcGravityForce(), PTCL::CalcGravityForce(), sph_system, dinfo); - #endif - PROBLEM::addExternalForce(sph_system, sysinfo); + if (PROBLEM::self_gravity) + grav_tree->calcForceAllAndWriteBack(PTCL::CalcGravityForce(), PTCL::CalcGravityForce(), sph_system, dinfo); + + PROBLEM::addExternalForce(sph_system, sysinfo); #pragma omp parallel for for(int i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i){ sph_system[i].finalKick(sysinfo.dt); diff --git a/src/mathfunc.h b/src/mathfunc.h index f7efe9f..dcf3bdf 100644 --- a/src/mathfunc.h +++ b/src/mathfunc.h @@ -1,9 +1,18 @@ #pragma once namespace math{ + // You probably created these functions, because you are concerned about speed in + // performance critical parts of your code? In this case you should replace the call + // to atan by just stating pi = 3.14159265358979323846. atan is actually an expensive + // function to call const PS::F64 pi = atan(1.0) * 4.0; - const PS::F64 NaN = + 0.0 / 0.0; - const PS::F64 VERY_LARGE_VALUE = 1.0e+30; + + // There is std::numeric_limits::signaling_NaN() for this purpose, or std::numeric_limits::quiet_NaN() + // if you do not want to throw an exception when this value is accessed + const PS::F64 NaN = std::numeric_limits::signaling_NaN(); + + // This is available as std::numeric_limits::max() + const PS::F64 VERY_LARGE_VALUE = std::numeric_limits::max(); template inline type plus(const type arg){ return (arg > 0) ? arg : 0; } diff --git a/src/param.h b/src/param.h index 9d65e9b..5d639a8 100644 --- a/src/param.h +++ b/src/param.h @@ -1,6 +1,11 @@ #pragma once +// It would be ideal to summarize all parameters in a single header file, so that a single file uniquely defines +// a certain model. Is it possible to move the content of this file into the init/GI.h file? It can then be duplicated and changed +// for every additional model setup + namespace PARAM{ + #ifdef PARTICLE_SIMULATOR_TWO_DIMENSION const short int Dim = 2; #else