Time averaging of calculated quantities

[makrandak]

Hello,

I am facing an interesting issue.

I am trying to compute long time averages of different terms that depend on fluctuating velocities and pressures. I implemented occa kernels for calculating these terms with the fluctuating quantities as input. I am trying to register these extra terms as fields in UDF_Setup(). Here is an example snippet.

  {
    auto mesh = nrs->meshV;

    auto avg = [&](const occa::memory& o_u)
    {
      return platform->linAlg->innerProd(mesh->Nlocal, mesh->o_Jw, o_u, platform->comm.mpiComm) / mesh->volume;
    };

    std::vector<dfloat> uAvg(mesh->dim);
    double pAvg;
    uAvg[0] = avg(nrs->o_U.slice(0 * nrs->fieldOffset));
    uAvg[1] = avg(nrs->o_U.slice(1 * nrs->fieldOffset));
    uAvg[2] = avg(nrs->o_U.slice(2 * nrs->fieldOffset));
    pAvg = avg(nrs->o_P.slice(0 * nrs->fieldOffset));

    auto o_uPrime = [&]()
    {
      poolDeviceMemory<dfloat> o_uPrime(nrs->NVfields * nrs->fieldOffset);
      o_uPrime.copyFrom(nrs->o_U, o_uPrime.size());
      for(int i = 0; i < mesh->dim; i++) {
          platform->linAlg->add(mesh->Nlocal, -uAvg[i], o_uPrime.slice(i * nrs->fieldOffset));
      }
      return o_uPrime;
    }();

    auto o_pPrime = [&]()
    {
      poolDeviceMemory<dfloat> o_pPrime(nrs->fieldOffset);
      o_pPrime.copyFrom(nrs->o_P, o_pPrime.size());
      platform->linAlg->add(mesh->Nlocal, -pAvg, o_pPrime);
      return o_pPrime;
    }();

    const auto epsilongVolumeAvg = [&]()
    {
      dfloat mue, rho;
      platform->options.getArgs("VISCOSITY", mue);
      platform->options.getArgs("DENSITY", rho);
      const auto nu = mue / rho;

      auto o_gradU = opSEM::strongGradVec(mesh, nrs->fieldOffset, o_uPrime);

      static poolDeviceMemory<dfloat> o_invAJwJw;
      if (!o_invAJwJw.isInitialized()) {
        o_invAJwJw.resize(mesh->Nlocal);
        platform->linAlg->axmyz(mesh->Nlocal, 1.0, mesh->o_invAJw, mesh->o_Jw, o_invAJwJw);
      }
      platform->linAlg->axmyMany(mesh->Nlocal, mesh->dim * mesh->dim, nrs->fieldOffset, 0, 1.0, mesh->o_invAJw, o_gradU);
      nrs->qqt->startFinish("+", o_gradU, nrs->fieldOffset);

      poolDeviceMemory<dfloat> o_eps(mesh->Nlocal);
      platform->linAlg->magSqrTensor(mesh->Nlocal, nrs->fieldOffset, o_gradU, o_eps);
      platform->linAlg->scale(mesh->Nlocal, nu, o_eps);
      return o_eps;
    }();

    const auto kVol = [&]()
    {
      poolDeviceMemory<dfloat> o_k(nrs->fieldOffset);
      //kEnergy(mesh->Nlocal, nrs->fieldOffset, o_uPrime, o_k);
      kEnergy(mesh->Nlocal, nrs->fieldOffset, nrs->o_U, o_k);
      return o_k;
    }();

    const auto dissipationTensor = [&]()
    {
      dfloat mue, rho;
      platform->options.getArgs("VISCOSITY", mue);
      platform->options.getArgs("DENSITY", rho);
      const auto nu = mue / rho;

      //auto o_gradU = opSEM::strongGradVec(mesh, nrs->fieldOffset, o_uPrime);
      auto o_gradU = opSEM::strongGradVec(mesh, nrs->fieldOffset, nrs->o_U);

      static poolDeviceMemory<dfloat> o_invAJwJw;
      if (!o_invAJwJw.isInitialized()) {
        o_invAJwJw.resize(mesh->Nlocal);
        platform->linAlg->axmyz(mesh->Nlocal, 1.0, mesh->o_invAJw, mesh->o_Jw, o_invAJwJw);
      }
      platform->linAlg->axmyMany(mesh->Nlocal, mesh->dim * mesh->dim, nrs->fieldOffset, 0, 1.0, mesh->o_invAJw, o_gradU);
      nrs->qqt->startFinish("+", o_gradU, nrs->fieldOffset);

      poolDeviceMemory<dfloat> o_disspTensor(9 * nrs->fieldOffset);
      o_disspTensor.copyFrom(o_gradU, o_disspTensor.size());
      dissipationTensorCalc(mesh->Nlocal, nrs->fieldOffset, nu, o_gradU, o_disspTensor);
      return o_disspTensor;
    }();

    std::vector<std::vector<deviceMemory<dfloat>>> tavgFields;

    deviceMemory<dfloat> o_u1(o_uPrime.slice(0 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_u2(o_uPrime.slice(1 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_u3(o_uPrime.slice(2 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_p(o_pPrime.slice(0 * nrs->fieldOffset , nrs->fieldOffset));

    /// Dissipation tensor
    deviceMemory<dfloat> o_dt11(dissipationTensor.slice(0 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_dt12(dissipationTensor.slice(1 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_dt13(dissipationTensor.slice(2 * nrs->fieldOffset , nrs->fieldOffset));

    deviceMemory<dfloat> o_dt21(dissipationTensor.slice(3 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_dt22(dissipationTensor.slice(4 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_dt23(dissipationTensor.slice(5 * nrs->fieldOffset , nrs->fieldOffset));

    deviceMemory<dfloat> o_dt31(dissipationTensor.slice(6 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_dt32(dissipationTensor.slice(7 * nrs->fieldOffset , nrs->fieldOffset));
    deviceMemory<dfloat> o_dt33(dissipationTensor.slice(8 * nrs->fieldOffset , nrs->fieldOffset));

    /// Register fluctuating velocities for time average
    tavgFields.push_back({o_u1});
    tavgFields.push_back({o_u2});
    tavgFields.push_back({o_u3});

    /// Register Reynolds stresses for time average
    tavgFields.push_back({o_u1, o_u1});
    tavgFields.push_back({o_u1, o_u2});
    tavgFields.push_back({o_u1, o_u3});
    tavgFields.push_back({o_u2, o_u1});
    tavgFields.push_back({o_u2, o_u2});
    tavgFields.push_back({o_u2, o_u3});
    tavgFields.push_back({o_u3, o_u1});
    tavgFields.push_back({o_u3, o_u2});
    tavgFields.push_back({o_u3, o_u3});

    /// Register dissipation tensor for time average
    tavgFields.push_back({o_dt11});
    tavgFields.push_back({o_dt12});
    tavgFields.push_back({o_dt13});
    tavgFields.push_back({o_dt21});
    tavgFields.push_back({o_dt22});
    tavgFields.push_back({o_dt23});
    tavgFields.push_back({o_dt31});
    tavgFields.push_back({o_dt32});
    tavgFields.push_back({o_dt33});

    /// Register rate of dissipation for time average
    tavgFields.push_back({o_epsilon});
    /// Register kinetic energy for time average
    tavgFields.push_back({o_kineticEnergy});

    tavg::setup(nrs->fieldOffset, tavgFields);
  }

So that you know, each term that is pushed back into tavgFields has a calculation associated with it. However, when I run this, all I get is zeros. My speculation is that these calculations go out of scope once UDF_Setup is done, and these calculations are never performed afterward (they need to be performed each timestep). Please correct me if I am incorrect in my understanding. If that is the case, is there a better way to register and compute time averages of such custom terms?

I really appreciate your help.

With sincere regards,
Makrand

[stgeke]

The variables won't go out of scope but I don't see how you update them at each time step as you declare them locally in UDF_Setup

[makrandak]

Yes, that is the problem. As they are defined in UDF_Setup(), they do not get update every timestep. Is there a way to reference them at each timestep and update them in UDF_Executestep()?

Right now, my hack solution is to define member variables in the nrs class and allocate them, and update them in UDF_Executestep(). Admittedly it is not a clean solution as it involves bypassing the udf design to define member variables in the nrs class. I would really appreciate your thoughts on a better way to do this.

With sincere regards,
Makrand

[stgeke]

Define dissipationTensor as a global variable in the udf.