Skip to content

Reply To: Implementing the Cumulant Collision Operator

Due to recent bot attacks we have changed the sign-up process. If you want to participate in our forum, first register on this website and then send a message via our contact form.

Forums OpenLB General Topics Implementing the Cumulant Collision Operator Reply To: Implementing the Cumulant Collision Operator

#10688
jmctighe
Participant

FBukreev,

I’ve written the following postprocessor for computing the full covariance matrix of velocity, making use of the already existing fields LATTICE_TIME and AVERAGE_VELOCITY:

struct TrackCovarVelocity {
  static constexpr OperatorScope scope = OperatorScope::PerCellWithParameters;

  int getPriority() const {
    return std::numeric_limits<int>::max();
  }
 
  using parameters = meta::list<descriptors::LATTICE_TIME>;

  template <typename CELL,typename PARAMETERS>
  void apply(CELL& cell, PARAMETERS& parameters) any_platform {
    using V = typename CELL::value_t;
    using DESCRIPTOR = typename CELL::descriptor_t;

    int iT = parameters.template get<descriptors::LATTICE_TIME>();
    auto uCovar = cell.template getField<COVARIANCE_VELOCITY>();
    if (iT > 1) { 
      auto uAvg = cell.template getField<AVERAGE_VELOCITY>();

      V u[DESCRIPTOR::d] {};
      cell.computeU(u);
      
      int k;
      for (int i=0; i < DESCRIPTOR::d; i++){
        for (int j = 0; j < DESCRIPTOR::d; j++){
          k = j + i * DESCRIPTOR::d;
         uCovar[k] = (uCovar[k]*(iT - 1) + (u[i] - uAvg[i]) * (u[j] - uAvg[j])*(iT/(iT-1))) / iT;
        }
      }     
    }
    else {
      for(int i=0; i<(DESCRIPTOR::d * DESCRIPTOR::d); i++){
        uCovar[i] = 0.;
      }
    }
    cell.template setField<COVARIANCE_VELOCITY>(uCovar);

  }
  
  
};

The postprocessor is working as expected. Thank you again for all your help!

Best,

jmctighe