Second declaration of divV should be removed in post::div_post

[dlansigan]

Description

divV is declared twice in post::div_post, at lines 394 and 416. Thus, when running a 3D simulation, the outer-scope divV is written, which is exactly 0.0. The fix would be to remove "double" in line 416.

Array<double> F;
double divV = 0.0;
double Jac = 0.0;
Array<double> ksix(nsd,nsd);

for (int g = 0; g < lM.nG; g++) {
  if (g == 0 || !lM.lShpF) {
    auto Nx_g = lM.Nx.slice(g);
    nn::gnn(eNoN, nsd, nsd, Nx_g, xl, Nx, Jac, ksix);
  }

  double w = lM.w(g)*Jac;
  N = lM.N.col(g);

  if ((cPhys == EquationType::phys_fluid) || (cPhys == EquationType::phys_CMM)) {
    Array<double> vx(nsd,nsd);

    if (nsd == 3) {
      for (int a = 0; a < eNoN; a++) {
        vx(0,0) = vx(0,0) + Nx(0,a)*yl(0,a);
        vx(1,1) = vx(1,1) + Nx(1,a)*yl(1,a);
        vx(2,2) = vx(2,2) + Nx(2,a)*yl(2,a);
      }
      double divV = vx(0,0) + vx(1,1) + vx(2,2);

    } else { 
      for (int a = 0; a < eNoN; a++) {
        vx(0,0) = vx(0,0) + Nx(0,a)*yl(0,a);
        vx(1,1) = vx(1,1) + Nx(1,a)*yl(1,a);
      }
      divV = vx(0,0) + vx(1,1);
    }

This does not occur in svFSI, so this is likely just a translation error. See the corresponding code in svFSI (divV is declared in line 1576):

IF ((cPhys.EQ.phys_fluid) .OR. (cPhys.EQ.phys_CMM)) THEN
   vx = 0._RKIND
   IF (nsd .EQ. 3) THEN
      DO a=1, eNoN
         vx(1,1) = vx(1,1) + Nx(1,a)*yl(1,a)
         vx(2,2) = vx(2,2) + Nx(2,a)*yl(2,a)
         vx(3,3) = vx(3,3) + Nx(3,a)*yl(3,a)
      END DO
      divV = vx(1,1) + vx(2,2) + vx(3,3)
   ELSE
      DO a=1, eNoN
         vx(1,1) = vx(1,1) + Nx(1,a)*yl(1,a)
         vx(2,2) = vx(2,2) + Nx(2,a)*yl(2,a)
      END DO
      divV = vx(1,1) + vx(2,2)
   END IF

Reproduction

Run any 3D fluid simulation (e.g., tests/cases/fluid/pipe_RCR_3d) and read dumped divergence values.

Expected behavior

Divergence should be close to zero but not exactly zero.

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[dlansigan]

Upon this fix, several of the test cases now give large divergence values. For example:

• fluid/driven_cavity_2d: min=-13.01, max=12.96
• fluid/iliac_artery: min=-186.92, max=170.80
• fluid/pipe_RCR_3d: min=-0.1523, max=0.1570

Are these magnitudes expected?

[aabrown100-git]

Thanks for reporting these. Could you also compute these normalized to the max velocity magnitude in the domain? And perhaps the RMS values over the domain?

I'm having trouble finding references for "normal" values of the velocity divergence in finite-element-based CFD. @dlansigan I think you can merge the pull request, and make a new issue about whether these divergence values are expected.