generate initial conditions forparticles with orientation

bfps/DNS.py

@@ -127,7 +127,7 @@ class DNS(_code):
                     template_class = '{0}<{1}>::'.format(self.dns_type, rnumber),
                     template_prefix = 'template '.format(rnumber),
                     just_declaration = True) + '\n\n')
-            if self.dns_type in ['NSVEparticles', 'NSVE_no_output', 'NSVEparticles_no_output']:
+            if self.dns_type in ['NSVEparticles', 'NSVE_no_output', 'NSVEparticles_no_output', 'NSVEparticlesP2P']:
                 outfile.write('template <typename rnumber> int NSVE<rnumber>::read_parameters(){return EXIT_SUCCESS;}\n')
                 outfile.write('template int NSVE<float>::read_parameters();\n')
                 outfile.write('template int NSVE<double>::read_parameters();\n\n')
@@ -377,7 +377,7 @@ class DNS(_code):
         assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0)
         assert (self.parameters['niter_todo'] % self.parameters['niter_out']  == 0)
         assert (self.parameters['niter_out']  % self.parameters['niter_stat'] == 0)
-        if self.dns_type in ['NSVEparticles_no_output', 'NSVEparticles']:
+        if self.dns_type in ['NSVEparticles_no_output', 'NSVEparticlesP2P', 'NSVEparticles']:
             assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
             assert (self.parameters['niter_out']  % self.parameters['niter_part'] == 0)
         _code.write_par(self, iter0 = iter0)
@@ -432,6 +432,9 @@ class DNS(_code):
             number_of_particles = 1
             for val in pbase_shape[1:]:
                 number_of_particles *= val
+        ncomponents = 3
+        if self.dns_type in ['NSVEparticlesP2P']:
+            ncomponents = 6
         with h5py.File(self.get_checkpoint_0_fname(), 'a') as ofile:
             s = 0
             ofile.create_group('tracers{0}'.format(s))
@@ -442,13 +445,13 @@ class DNS(_code):
                     shape = (
                         (self.parameters['tracers{0}_integration_steps'.format(s)],) +
                         pbase_shape +
-                        (3,)),
+                        (ncomponents,)),
                     dtype = np.float)
             ofile['tracers{0}/state'.format(s)].create_dataset(
                     '0',
                     shape = (
                         pbase_shape +
-                        (3,)),
+                        (ncomponents,)),
                     dtype = np.float)
         return None
     def job_parser_arguments(
@@ -621,6 +624,17 @@ class DNS(_code):
         self.parameters_to_parser_arguments(
                 parser_NSVEp2,
                 self.NSVEp_extra_parameters)
+
+        parser_NSVEp2p = subparsers.add_parser(
+                'NSVEparticlesP2P',
+                help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
+        self.simulation_parser_arguments(parser_NSVEp2p)
+        self.job_parser_arguments(parser_NSVEp2p)
+        self.particle_parser_arguments(parser_NSVEp2p)
+        self.parameters_to_parser_arguments(parser_NSVEp2p)
+        self.parameters_to_parser_arguments(
+                parser_NSVEp2p,
+                self.NSVEp_extra_parameters)
         return None
     def prepare_launch(
             self,
@@ -656,7 +670,7 @@ class DNS(_code):
         self.dns_type = opt.DNS_class
         self.name = self.dns_type + '-' + self.fluid_precision + '-v' + bfps.__version__
         # merge parameters if needed
-        if self.dns_type in ['NSVEparticles', 'NSVEparticles_no_output']:
+        if self.dns_type in ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']:
             for k in self.NSVEp_extra_parameters.keys():
                 self.parameters[k] = self.NSVEp_extra_parameters[k]
         if type(extra_parameters) != type(None):
@@ -690,7 +704,7 @@ class DNS(_code):
             # hardcoded FFTW complex representation size
             field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize
             checkpoint_size = field_size
-            if self.dns_type in ['NSVEparticles', 'NSVEparticles_no_output']:
+            if self.dns_type in ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']:
                 rhs_size = self.parameters['tracers0_integration_steps']
                 if type(opt.tracers0_integration_steps) != type(None):
                     rhs_size = opt.tracers0_integration_steps
@@ -726,14 +740,24 @@ class DNS(_code):
             nn = self.parameters['nparticles']
             cc = int(0)
             batch_size = int(1e6)
+            def get_random_phases(npoints):
+                return np.random.random(
+                            (npoints, 3))*2*np.pi
+            def get_random_versors(npoints):
+                bla = np.random.normal(
+                        size = (npoints, 3))
+                bla  /= np.sum(bla**2, axis = 1)[:, None]
+                return bla
             while nn > 0:
                 if nn > batch_size:
-                    dset[cc*batch_size:(cc+1)*batch_size] = np.random.random(
-                            (batch_size, 3))*2*np.pi
+                    dset[cc*batch_size:(cc+1)*batch_size, :3] = get_random_phases(batch_size)
+                    if dset.shape[1] == 6:
+                        dset[cc*batch_size:(cc+1)*batch_size, 3:] = get_random_versors(batch_size)
                     nn -= batch_size
                 else:
-                    dset[cc*batch_size:cc*batch_size+nn] = np.random.random(
-                            (nn, 3))*2*np.pi
+                    dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn)
+                    if dset.shape[1] == 6:
+                        dset[cc*batch_size:cc*batch_size+nn, 3:] = get_random_versors(nn)
                     nn = 0
                 cc += 1
         return None
@@ -943,7 +967,7 @@ class DNS(_code):
             #            particle_initial_condition[..., 2] += onedarray[None, :, None, None]
             self.write_par(
                     particle_ic = None)
-            if self.dns_type in ['NSVEparticles', 'NSVEparticles_no_output']:
+            if self.dns_type in ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']:
                 self.generate_particle_data(opt = opt)
         self.run(
                 nb_processes = opt.nb_processes,