[do not merge] Make some higher level comments about the code

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

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 Ptcl> class Problem{
 	Problem(){
 	}

src/init/GI.h

@@ -1,11 +1,12 @@
-#define SELF_GRAVITY
 #define FLAG_GI
 #ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
 #error
 #endif
 template <class Ptcl> class GI : public Problem<Ptcl>{
 	public:
 	static const double END_TIME = 1.0e+4;
+	static const bool self_gravity = true;
+
 	static void setupIC(PS::ParticleSystem<Ptcl>& 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

src/main.cpp

@@ -5,6 +5,7 @@
 #include <cstdlib>
 #include <iostream>
 #include <vector>
+#include <memory>
 
 #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<PTCL::RealPtcl> PROBLEM;
 	//////////////////
 	//Create vars.
@@ -56,27 +68,37 @@ int main(int argc, char* argv[]){
 	PS::TreeForForceShort<PTCL::RESULT::Dens , PTCL::EPI::Dens , PTCL::EPJ::Dens >::Gather   dens_tree;
 	PS::TreeForForceShort<PTCL::RESULT::Drvt , PTCL::EPI::Drvt , PTCL::EPJ::Drvt >::Gather   drvt_tree;
 	PS::TreeForForceShort<PTCL::RESULT::Hydro, PTCL::EPI::Hydro, PTCL::EPJ::Hydro>::Symmetry hydr_tree;
-	#ifdef SELF_GRAVITY
-	PS::TreeForForceLong <PTCL::RESULT::Grav , PTCL::EPI::Grav , PTCL::EPJ::Grav >::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<PS::TreeForForceLong <PTCL::RESULT::Grav , PTCL::EPI::Grav , PTCL::EPJ::Grav >::Monopole> grav_tree;
+	if (PROBLEM::self_gravity)
+	  grav_tree.reset(new PS::TreeForForceLong <PTCL::RESULT::Grav , PTCL::EPI::Grav , PTCL::EPJ::Grav >::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::EPJ::Grav>(), PTCL::CalcGravityForce<PS::SPJMonopole>(), sph_system, dinfo);
-	#endif
-	sysinfo.dt = getTimeStepGlobal<PTCL::RealPtcl>(sph_system);
+        if (PROBLEM::self_gravity)
+          grav_tree->calcForceAllAndWriteBack(PTCL::CalcGravityForce<PTCL::EPJ::Grav>(), PTCL::CalcGravityForce<PS::SPJMonopole>(), sph_system, dinfo);
+
+        sysinfo.dt = getTimeStepGlobal<PTCL::RealPtcl>(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::EPJ::Grav>(), PTCL::CalcGravityForce<PS::SPJMonopole>(), sph_system, dinfo);
-		#endif
-		PROBLEM::addExternalForce(sph_system, sysinfo);
+	        if (PROBLEM::self_gravity)
+	          grav_tree->calcForceAllAndWriteBack(PTCL::CalcGravityForce<PTCL::EPJ::Grav>(), PTCL::CalcGravityForce<PS::SPJMonopole>(), 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);

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<PS::F64>::signaling_NaN() for this purpose, or std::numeric_limits<PS::F64>::quiet_NaN()
+	// if you do not want to throw an exception when this value is accessed
+	const PS::F64 NaN = std::numeric_limits<PS::F64>::signaling_NaN();
+
+	// This is available as std::numeric_limits<PS::F64>::max()
+	const PS::F64 VERY_LARGE_VALUE = std::numeric_limits<PS::F64>::max();
 	template <typename type> inline type plus(const type arg){
 		return (arg > 0) ? arg : 0;
 	}

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