broken. working on tracers

bfps/test.py

@@ -38,12 +38,14 @@ class convergence_test(bfps.code):
         self.includes += '#include "tracers.hpp"\n'
         self.variables += ('double t;\n' +
                            'FILE *stat_file;\n'
+                           'FILE *traj_file;\n'
                            'double stats[2];\n')
         self.variables += self.cdef_pars()
         self.definitions += self.cread_pars()
         self.definitions += """
                 //begincpp
-                void do_stats(fluid_solver<float> *fsolver)
+                void do_stats(fluid_solver<float> *fsolver,
+                              tracer<float> *tracer)
                 {
                     fsolver->compute_velocity(fsolver->cvorticity);
                     stats[0] = .5*fsolver->correl_vec(fsolver->cvelocity,  fsolver->cvelocity);
@@ -53,6 +55,7 @@ class convergence_test(bfps.code):
                         fwrite((void*)&fsolver->iteration, sizeof(int), 1, stat_file);
                         fwrite((void*)&t, sizeof(double), 1, stat_file);
                         fwrite((void*)stats, sizeof(double), 2, stat_file);
+                        fwrite((void*)tracer->state, sizeof(double), tracer->array_size);
                     }
                 }
                 //endcpp
@@ -78,6 +81,8 @@ class convergence_test(bfps.code):
                     {
                         sprintf(fname, "%s_stats.bin", simname);
                         stat_file = fopen(fname, "wb");
+                        sprintf(fname, "%s_traj.bin", ps->name);
+                        traj_file = fopen(fname, "wb");
                     }
                 fs->read('v', 'c');
                 fs->low_pass_Fourier(fs->cvorticity, 3, fs->kM);
@@ -87,20 +92,25 @@ class convergence_test(bfps.code):
                         simname, fs,
                         32, 2,
                         fs->ru);
+                ps->dt = dt;
                 fs->compute_velocity(fs->cvorticity);
                 fftwf_execute(*((fftwf_plan*)fs->c2r_velocity));
-                ps->transfer_data();
+                ps->update_field();
                 t = 0.0;
-                do_stats(fs);
+                do_stats(fs, ps);
                 fs->write('u', 'r');
                 fs->write('v', 'r');
                 for (; fs->iteration < iter0 + niter_todo;)
                 {
                     fs->step(dt);
                     t += dt;
-                    do_stats(fs);
+                    ps->update_field();
+                    ps->Euler();
+                    ps->synchronize();
+                    do_stats(fs, ps);
                 }
                 fclose(stat_file);
+                fclose(traj_file);
                 fs->write('v', 'c');
                 fs->write('v', 'r');
                 fs->write('u', 'r');

src/slab_field_particles.cpp

@@ -214,6 +214,20 @@ void slab_field_particles<rnumber>::rFFTW_to_buffered(rnumber *src, rnumber *dst
             mpirequest);
     delete mpirequest;
 }
+
+
+
+template <class rnumber>
+void slab_field_particles<rnumber>::Euler()
+{
+    double *y = fftw_alloc_real(this->array_size);
+    this->get_rhs(this->state, y);
+    for (int p=0; p<this->nparticles; p++)
+        for (int i=0; i<this->ncomponents; i++)
+            this->state[p*this->ncomponents+i] += this->dt*y[p*this->ncomponents+i];
+    fftw_free(y);
+}
+
 /*****************************************************************************/
 /* finally, force generation of code for single precision                    */
 template class slab_field_particles<float>;

src/slab_field_particles.hpp

@@ -84,6 +84,9 @@ class slab_field_particles
         void synchronize();
         void rFFTW_to_buffered(rnumber *src, rnumber *dst);
         ptrdiff_t buffered_local_size();
+
+        /* solvers */
+        void Euler();
 };
 
 

src/tracers.cpp

@@ -32,11 +32,12 @@ void tracers<rnumber>::jump_estimate(double *jump)
 template <class rnumber>
 void tracers<rnumber>::get_rhs(double *x, double *y)
 {
-    std::fill_n(y, this->array_size, 0.0);
+    /* get grid coordinates */
+    /* perform interpolation */
 }
 
 template<class rnumber>
-void tracers<rnumber>::transfer_data(bool clip_on)
+void tracers<rnumber>::update_field(bool clip_on)
 {
     if (clip_on)
         clip_zero_padding(this->fs->rd, this->source_data, 3);

src/tracers.hpp

@@ -44,7 +44,7 @@ class tracers:public slab_field_particles<rnumber>
                 rnumber *SOURCE_DATA);
         ~tracers();
 
-        void transfer_data(bool clip_on = true);
+        void update_field(bool clip_on = true);
         virtual void get_rhs(double *x, double *rhs);
         virtual void jump_estimate(double *jump_length);
 };