Fix bug of adding particles in the middle of simulations for PeTar (#988)

Author: lwang-astro

src/amuse/community/petar/Makefile

@@ -40,7 +40,7 @@ LDFLAGS  += -lm $(MUSE_LD_FLAGS)
 
 OBJS = interface.o 
 
-FILELIST=interface.cc interface.py interface.h
+DIRLIST=src/PeTar src/SDAR src/FDPS
 
 CODELIB = 
 
@@ -60,18 +60,16 @@ distclean: clean
 
 DOWNLOAD_FROM_WEB = $(PYTHON) ./download.py
 
-download:
-	$(RM) -Rf .pc
-	$(RM) -Rf src
-	mkdir src
+$(DIRLIST):
+#	$(RM) -Rf .pc
+#	$(RM) -Rf src
+#       mkdir src
 	$(DOWNLOAD_FROM_WEB)
 
 __init__.py:
 	touch $@
 
-$(FILELIST): |download
-
-src/PeTar/src/get_version.hpp: src/PeTar/src/get_version.hpp.in |src/PeTar
+src/PeTar/src/get_version.hpp: |src/PeTar
 	sed 's/@VERSION@/'`cat src/PeTar/VERSION`'_'`cat src/SDAR/VERSION`'/g' src/PeTar/src/get_version.hpp.in >src/PeTar/src/get_version.hpp
 
 test_interface: $(OBJS) test_interface.o
@@ -86,5 +84,5 @@ worker_code.h: interface.py
 petar_worker: worker_code.cc worker_code.h $(CODELIB) $(OBJS) interface.h
 	$(MPICXX) $(CXXFLAGS) $(SC_FLAGS) $(AM_CFLAGS) $(LDFLAGS) $< $(OBJS) $(CODELIB) -o $@ $(SC_MPI_CLIBS) $(LIBS) $(AM_LIBS)
 
-interface.o: interface.cc src/PeTar/src/get_version.hpp 
+interface.o: interface.cc src/PeTar/src/get_version.hpp |$(DIRLIST)
 	$(MPICXX) $(CXXFLAGS) $(SC_FLAGS) -c -o $@ $< 

src/amuse/community/petar/download.py

@@ -47,21 +47,25 @@ def unpack_downloaded_file(self, filename, name, version):
         print("done")
 
     def start(self):
-        if os.path.exists('src'):
-            counter = 0
-            while os.path.exists('src.{0}'.format(counter)):
-                counter += 1
-                if counter > 100:
-                    print("too many backup directories")
-                    break
-            os.rename('src', 'src.{0}'.format(counter))
-
-        os.mkdir('src')
+        if not os.path.exists('src'):
+            os.mkdir('src')
+#            if not update_flag:
+#                return
+#            counter = 0
+#            while os.path.exists('src.{0}'.format(counter)):
+#                counter += 1
+#                if counter > 100:
+#                    print("too many backup directories")
+#                    break
+#            os.rename('src', 'src.{0}'.format(counter))
 
         for i, url_template in enumerate(self.url_template):
             url = url_template.format(version=self.version[i])
             filename = self.filename_template.format(version=self.version[i])
             filepath = os.path.join(self.src_directory(), filename)
+            if os.path.exists('src/'+self.name[i]):
+                print("src/%s exist" % self.name[i])
+                continue
             print(
                 "downloading version", self.version[i],
                 "from", url, "to", filename

src/amuse/community/petar/interface.cc

@@ -23,6 +23,7 @@ extern "C" {
 
     // flags
     static bool particle_list_change_flag=true;
+    static bool initial_particle_flag=false;
 
     // common
 
@@ -102,6 +103,8 @@ extern "C" {
         if(ptr->my_rank==0) std::cout<<"PETAR: recommit_parameters start\n";
 #endif
         ptr->input_parameters.n_glb.value = ptr->stat.n_real_glb;
+        ptr->input_parameters.update_changeover_flag = true;
+        ptr->input_parameters.update_rsearch_flag = true;
         ptr->initialParameters();
         ptr->initial_step_flag = false;
 
@@ -123,6 +126,15 @@ extern "C" {
 
     int new_particle(int* index_of_the_particle,  double mass, double x, double y, double z, double vx, double vy, double vz, double radius) {
         // if not yet initial the system
+#ifdef INTERFACE_DEBUG_PRINT
+        if (!ptr->read_data_flag && ptr->my_rank==0) {
+            std::cout<<"New particle, rank "<<ptr->my_rank<<std::endl;
+            std::cout<<std::setw(20)<<"ID";
+            ParticleBase::printColumnTitle(std::cout);
+            std::cout<<std::endl;
+        }
+#endif
+
         ptr->read_data_flag = true;
 
         PS::S64 id_offset = ptr->input_parameters.id_offset.value;
@@ -152,10 +164,25 @@ extern "C" {
             p.rank_org = ptr->my_rank;
             p.adr = n_loc;
 
+            if (initial_particle_flag) {
+                PS::F64 m_fac = p.mass*Ptcl::mean_mass_inv;
+                PS::F64 r_out = ptr->input_parameters.r_out.value;
+                PS::F64 r_in = r_out*ptr->input_parameters.ratio_r_cut.value;
+                p.changeover.setR(m_fac, r_in, r_out);
+                p.calcRSearch(ptr->input_parameters.dt_soft.value);
+            }
+            
             ptr->system_soft.addOneParticle(p);
 
             ptr->stat.n_real_loc++;
             *index_of_the_particle = p.id;
+#ifdef INTERFACE_DEBUG_PRINT
+            std::cout<<std::setprecision(14);
+            std::cout<<std::setw(20)<<p.id;
+            p.ParticleBase::printColumn(std::cout);
+            std::cout<<std::endl;
+#endif
+
         }
         ptr->stat.n_real_glb++;
 
@@ -548,7 +575,7 @@ extern "C" {
 
     int evolve_model(double time_next) {
 #ifdef INTERFACE_DEBUG_PRINT
-        if(ptr->my_rank==0) std::cout<<"PETAR: evolve models start\n";
+        if(ptr->my_rank==0) std::cout<<"PETAR: evolve models to "<<time_next<< "start\n";
 #endif
 
         if (ptr->stat.n_real_glb==0) {// escape if no particle
@@ -731,6 +758,7 @@ extern "C" {
         ptr->initialStep();
         ptr->reconstructIdAdrMap();
         particle_list_change_flag = false;
+        initial_particle_flag = true;
 #ifdef INTERFACE_DEBUG_PRINT
         if(ptr->my_rank==0) std::cout<<"PETAR: commit_particles end\n";
 #endif
@@ -859,6 +887,17 @@ extern "C" {
         return 0;
     }
 
+    int set_output_step(double dt_snap) {
+        ptr->input_parameters.dt_snap.value = dt_snap;
+        return 0;
+    }
+
+    int get_output_step(double* dt_snap) {
+        *dt_snap = ptr->input_parameters.dt_snap.value;
+        return 0;
+    }
+
+
     //// set gravitational constant
     //int set_gravitational_constant(double G) {
     //    ptr->input_parameters.unit_set.value=-1;

src/amuse/community/petar/interface.h

@@ -87,6 +87,10 @@ int get_tree_step(double * dt_soft);
 
 int set_tree_step(double dt_soft);
 
+int get_output_step(double * dt_output);
+
+int set_output_step(double dt_output);
+
 int get_kinetic_energy(double * kinetic_energy);
 
 int get_potential_energy(double * potential_energy);

src/amuse/community/petar/interface.py

@@ -356,6 +356,48 @@ def set_tree_step():
         """
         return function
 
+    @legacy_function
+    def get_output_step():
+        """
+        Get dt_output, the PeTar internal output time step (0.125)
+        """
+        function = LegacyFunctionSpecification()
+        function.addParameter(
+            'dt_output', dtype='float64', direction=function.OUT,
+            description=(
+                "dt_output, the PeTar internal output time step (0.125)"
+            )
+        )
+        function.result_type = 'int32'
+        function.result_doc = """
+        0 - OK
+            the parameter was retrieved
+        -1 - ERROR
+            could not retrieve parameter
+        """
+        return function
+
+    @legacy_function
+    def set_output_step():
+        """
+        Set dt_output, the PeTar internal output time step (0.125)
+        """
+        function = LegacyFunctionSpecification()
+        function.addParameter(
+            'dt_output', dtype='float64', direction=function.IN,
+            description=(
+                "dt_output, the PeTar internal output time step (0.125)"
+            )
+        )
+        function.result_type = 'int32'
+        function.result_doc = """
+        0 - OK
+            the parameter was set
+        -1 - ERROR
+            could not set parameter
+        """
+        return function
+
 
 class petar(GravitationalDynamics, GravityFieldCode):
 
@@ -443,6 +485,14 @@ def define_parameters(self, handler):
             default_value=0.0 | nbody_system.time
         )
 
+        handler.add_method_parameter(
+            "get_output_step",
+            "set_output_step",
+            "dt_output",
+            "PeTar internal output time step (0.125)",
+            default_value=0.125 | nbody_system.time
+        )
+
     def define_methods(self, handler):
         GravitationalDynamics.define_methods(self, handler)
         self.stopping_conditions.define_methods(handler)
@@ -542,6 +592,25 @@ def define_methods(self, handler):
             )
         )
 
+        handler.add_method(
+            "set_output_step",
+            (
+                nbody_system.time,
+            ),
+            (
+                handler.ERROR_CODE,
+            )
+        )
+
+        handler.add_method(
+            "get_output_step",
+            (),
+            (
+                nbody_system.time,
+                handler.ERROR_CODE,
+            )
+        )
+
     def define_particle_sets(self, handler):
         GravitationalDynamics.define_particle_sets(self, handler)
         self.stopping_conditions.define_particle_set(handler)