From ddc5c789d6db3210ce612784c78164e333ec3b79 Mon Sep 17 00:00:00 2001
From: evan1022 <evan.scannapieco@asu.edu>
Date: Fri, 8 Mar 2024 12:53:44 -0800
Subject: [PATCH] Update turbulence.cpp

---
 src/pgen/turbulence.cpp | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp
index 23f101ee..f8138f39 100644
--- a/src/pgen/turbulence.cpp
+++ b/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