Update turbulence.cpp

src/pgen/turbulence.cpp

@@ -536,5 +536,34 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
   auto state_dist = few_modes_ft.GetDistState();
   pin->SetString("problem/turbulence", "state_dist", state_dist);
 }
+\\
+\\ COMPUTE THE AVERAGE OVER PRODUCTS OF S VALUES
+\\ PHILIPP: PLEASE FIX THIS
+\\
+auto correlation_s = tracers_pkg->Param<std::vector<Real>>("correlation_s");
+auto correlation_sdot = tracers_pkg->Param<std::vector<Real>>("correlation_sdot");
+
+Kokkos::parallel_reduce(
+      "",
+      Kokkos::MDRangePolicy<Kokkos::Rank<1>>(n),
+      KOKKOS_LAMBDA(const int n, Real &corr_s, Real &corr_sdot) {
+        auto tracers_pkg = pmb->packages.Get("tracers");
+        const auto n_lookback = tracers_pkg->Param<int>("n_lookback");
+        auto idx = n_lookback - 1;
+        auto &s = swarm->Get<Real>("s").Get();
+        auto &sdot = swarm->Get<Real>("sdot").Get();
+        auto corr_s;
+        auto corr_sdot;
+
+        for (i = 0; i < idx; i++) {
+          corr_s[i]   =s(0, n)   *s(i,n);  \\PHILIPP: HOW DO YOU DO THIS INDEXING?
+          corr_sdot[i]=sdot(0, n)*sdot(i,n);
+        }      
+      },
+      Kokkos::Add<Real>(correlation_s,correlation_sdot));
+      correlation_s=correlation_s/N;
+      correlation_sdot=correlation_sdot/N;
+
+
 
 } // namespace turbulence