lwang-astro / PeTar

PeTar is a high-performance N-body code for modelling the evolution of star clusters and tidal streams, including the effect of galactic potential, dynamics of binary and hierarchical system, single and binary stellar evolution.
MIT License
71 stars 19 forks source link

Makefile for AMUSE should not hardcode architecture #28

Closed rieder closed 2 years ago

rieder commented 2 years ago

Currently the Makefile for the AMUSE version of PeTar is hardcoded for x86:

# arch
CXXFLAGS += -march=core-avx2
CXXFLAGS += -D INTRINSIC_X86
CXXFLAGS += -D USE_SIMD
CXXFLAGS += -D DIV_FIX

This prevents it from building on other architectures (e.g. arm64).

However, simply commenting these out doesn't seem to work, as compiling then runs into errors:

interface.cc:17:12: error: 'CalcForcePPSimd' does not name a type
   17 |     static CalcForcePPSimd<ParticleBase,FPSoft> fcalc; // Force calculator
      |            ^~~~~~~~~~~~~~~
lwang-astro commented 2 years ago

I see. The force calculation between particle and particle for bridge does not include a no simd version. Now I have added some update but not yet checked whether it works since my current working computer has not yet installed amuse. Could you try to implement the following change and see whether it provides the correct result for the bridge force? My test code in amuse-interface can also check this. By the way, I found that for the bridge force, I fix the softening length to zero. Is it fine?

diff --git a/amuse-interface/interface.cc b/amuse-interface/interface.cc
index 3c96fdf..d4bd8a0 100644
--- a/amuse-interface/interface.cc
+++ b/amuse-interface/interface.cc
@@ -14,7 +14,11 @@ extern "C" {

     static PeTar* ptr=NULL;
     static double time_start = 0.0;
+#ifdef USE_SIMD
     static CalcForcePPSimd<ParticleBase,FPSoft> fcalc; // Force calculator
+#else
+    static CalcForcePPNoSimd<ParticleBase,FPSoft> fcalc;
+#endif
     static int n_particle_in_interrupt_connected_cluster_glb; //

     // flags
diff --git a/src/soft_force.hpp b/src/soft_force.hpp
index 7935e65..c97be3a 100644
--- a/src/soft_force.hpp
+++ b/src/soft_force.hpp
@@ -199,6 +199,42 @@ struct CalcForceEpSpQuadNoSimd{
     }
 };

+template<class Tpi, class Tpj>
+struct CalcForcePPNoSimd {
+    void operator () (const Tpi * ep_i,
+                      const PS::S32 n_ip,
+                      const Tpj * ep_j,
+                      const PS::S32 n_jp,
+                      ForceSoft * force){
+      //const PS::F64 eps2 = EPISoft::eps * EPISoft::eps;
+        const PS::F64 eps2 = 0;
+        const PS::F64 G = ForceSoft::grav_const;
+        for(PS::S32 i=0; i<n_ip; i++){
+            PS::F64vec xi = ep_i[i].pos;
+            PS::F64vec ai = 0.0;
+            PS::F64 poti = 0.0;
+            for(PS::S32 j=0; j<n_jp; j++){
+                PS::F64vec rij = xi - ep_j[j].getPos();
+                PS::F64 r3_inv = rij * rij + eps2;
+                PS::F64 r_inv = 1.0/sqrt(r3_inv);
+                r3_inv = r_inv * r_inv;
+                r_inv *= ep_j[j].getCharge();
+                r3_inv *= r_inv;
+                ai -= r3_inv * rij;
+                poti -= r_inv;
+            }
+            force[i].acc += G*ai;
+            force[i].pot += G*poti;
+#ifdef NAN_CHECK_DEBUG
+            assert(!std::isnan(ai[0]));
+            assert(!std::isnan(ai[1]));
+            assert(!std::isnan(ai[2]));
+            assert(!std::isnan(poti));
+#endif
+        }
+    }
+};
+
 #ifdef USE_SIMD
rieder commented 2 years ago

It's not working yet, I get this error:

src/PeTar/src/soft_force.hpp:221:34: error: 'const class FPSoft' has no member named 'getCharge'
  221 |                 r_inv *= ep_j[j].getCharge();
      |                          ~~~~~~~~^~~~~~~~~
rieder commented 2 years ago

Bridge in AMUSE can work fine with zero softening.

rieder commented 2 years ago

If it's of any use I would have time next week or so to meet and discuss this / help with installing AMUSE if needed.

lwang-astro commented 2 years ago

I think these two corrections for soft_force.hpp:CalcForcePPNoSimd can fix the issue.

                PS::F64vec rij = xi - ep_j[j].pos;
                r_inv *= ep_j[j].mass;
rieder commented 2 years ago

yes now it compiles. I will test later.

rieder commented 2 years ago

test_interface.py runs fine!