ComputeCurl routines introduces errors in pressure at partition boundaries

[shaering]

Communications portion of ComputeCurl routines in utils causes partition imprinting on pressure solution. Gradients appear to be low at partition boundaries. Removing of communication portion

  // Count the zones globally.                                                                                                                                                  
  GroupCommunicator &gcomm = u.ParFESpace()->GroupComm();                                                                                                                       
  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);                                                                                                                    
  gcomm.Bcast(zones_per_vdof);                                                                                                                                                  

  // Accumulate for all vdofs.                                                                                                                                                  
  gcomm.Reduce<double>(cu.GetData(), GroupCommunicator::Sum);                                                                                                                   
  gcomm.Bcast<double>(cu.GetData()); 

fixes the issue. Not apparent as to why these sums introduce error into following mean calculations.

[trevilo]

I was unable to reproduce any partition imprinting on the pressure. After discussion with @shaering , it became clear that the partition imprinting was first observed while implementing the dissipation calculation. I was able to observe a partition imprinting in the dissipation when using the dev-dissipationStat branch with a variant of the gradient calculation. Specifically, switching from

      // gradient of u'
      setScalarFromVector(tmpR1_, 0, &tmpR0_);
      G_op_->Mult(tmpR0_, gradU_);
      setScalarFromVector(tmpR1_, 1, &tmpR0_);
      G_op_->Mult(tmpR0_, gradV_);
      if (dim_ == 3) {
        setScalarFromVector(tmpR1_, 2, &tmpR0_);
        G_op_->Mult(tmpR0_, gradW_);
      }

      Mv_inv_->Mult(gradU_, tmpR1_);
      gradU_ = tmpR1_;
      Mv_inv_->Mult(gradV_, tmpR1_);
      gradV_ = tmpR1_;
      if (dim_ == 3) {
        Mv_inv_->Mult(gradW_, tmpR1_);
        gradW_ = tmpR1_;
      }

      gradU_gf_->SetFromTrueDofs(gradU_);
      gradV_gf_->SetFromTrueDofs(gradV_);
      if (dim_ == 3) gradW_gf_->SetFromTrueDofs(gradW_);

to

      Up_gf_->SetFromTrueDofs(tmpR1_);
      vectorGrad3D(*epsi_gf_, *Up_gf_, *gradU_gf_, *gradV_gf_, *gradW_gf_);

in Tomboulides::computeDissipation shows clear evidence of a problem at partition boundaries in the resulting dissipation. See the top plot below, which shows the dissipation computed from a low Mach channel flow test case.

However, this problem results from an apparent bug in the scalarGrad3D function that is called from vectorGrad3D. Specifically, in scalarGrad3D, the GroupCommunicator object from the scalar input function is used to communicate the the gradient vector (see the code starting here). If we replace this GroupCommunicator with one obtained from the gradient ParGridFunction, then there is no evidence of the partitioning in the result (as shown in the bottom plot below).

This problem---i.e., a mismatch between the GroupCommunicator and the data being communicated---should not affect the curl, since the velocity vector and the curl approximation are using the same finite element space. In addition to the pressure, the curl of the velocity does not show any obvious evidence of the partitioning. Thus, as far as I can tell at the moment, the problem is restricted to scalarGrad3D, which is not being used in the code at this point, at least on main, as far as I know, but should nonetheless be fixed (or maybe just eliminated unless there is a compelling reason to keep it).