Getting the coordinates of quadrature points for the CHT case

[joneuhauser]

tl;dr: The coordinates I'm getting only make sense for a single MPI rank, and I can't figure out why.

As an example, I'll use the conj_ht example.

I've added the following to UDF_ExecuteStep to debug the coordinates:

  auto mesh = nrs->cds->mesh[0];
  std::ofstream outfile(std::to_string(platform->comm.mpiCommSize) + "coordinates" + std::to_string(platform->comm.mpiRank) + ".txt", std::ios::trunc);

  std::cout << mesh->Nelements << std::endl;
  for (int id = 0; id < mesh->Nelements * mesh->Np; id++)
  {
    outfile << id << ";" << mesh->x[id] << ";" << mesh->y[id] << ";" << mesh->z[id] << std::endl;
  }
  exit(0);

In the case of conj_ht, there are 96 (nx = 8, ny = 4, nz =3) elements in the meshV and 96 (nx = 8, ny = 2, nz = 3, times 2) elements in the solid domain, i.e. 192 elements in the nrs->cds->mesh[0]. The fluid has y coordinates between 0 and 1, and the solid has y coordinates <0 and >1.

If I look at the resulting coordinates files:

• For nrsmpi conj_ht 1, there are mesh->Nelements * mesh->Np entries in the file (good), and the second half of the file has coordinates inside the solid (good).
• For nrsmpi conj_ht 2, there are mesh->Nelements * mesh->Np / 2 entries in each file (good), but there are no coordinates inside the solid anywhere (i.e. not a single minus sign) (bad).

How do I get the coordinates inside the solid for a CHT case?

Context: I'm implementing averaging in the homogeneous direction and wanted to validate my node mapping.

[stgeke]

Did you try to outfile.close() before calling exit(0)?

[stgeke]

I guess that's not the issue here as each file has the correct number of entries.

[stgeke]

Could it be that rank0/rank1 owns all the solid element top/buttom. This would explain why one file doesn't contain any negative coordinates.

[joneuhauser]

Thanks for your answer! Neither of the two created files (for 2 ranks) contains any negative coordinates. All ranks have a part of the fluid domain, and they all show up in the files, and the global element id nek::lglel is as expected. Not so much for the solid domain

If I join up all of the generated files - they should contain the coordinates of all nodes in the domain, at least that's what I think my code does - and plot a histogram of y coordinates, here's what I get:

The line for ranks=1 looks as expected (double the amount of nodes at all of the element interfaces). For ranks=2 and ranks=3, no negative coordinates (and coordinates > 1) show up. The integral under all histograms is identical (i.e. the same total number of nodes is printed).

I'm working with an unmodified v22.0.0 (8ee9c38) on MVAPICH 2.3.7 (CUDA support) with the conj_ht example - only modification is displayed above (and I've reduced the polynomial order to 2 for easier debugging), and the outcome is identical for OCCA mode Cuda or Serial.

[stgeke]

Can you dump just the solid elements in the t-mesh:

  for (int e = 0; e < mesh->Nelements; e++)
  {
    if(mesh->elementInfo[e] == 1) { // solid
      for (int i = 0; i < mesh->Np; i++) {
        const auto id = i + e*mesh->Np;
        outfile << mesh->x[id] << "  " << mesh->y[id] << "  " << mesh->z[id] << std::endl;
      }
    }
  }

I just tested it (on 2 ranks) using the conj_ht example and it works fine (outcome is as expected).
btw: I used latest next not v22.0 but this shouldn't matter

[joneuhauser]

Yes, this is with an otherwise unmodified conj_ht, I even recompiled my MPI and NekRS just to be sure. Alright, I'll try the next branch!

[stgeke]

The y-coord of the solid elements in the conj_ht range from y = -0.5 .. 0 (buttom) and 1 .. 1.5 (top)

[joneuhauser]

Yes. So there should be negative y coordinates in one of those files (as there are for the ranks=1 case).

With 23.0.3 all seems well, there are negative y coordinates. Here's the same plot as above for 23.0.3:

(curves coincide, so all is well).

I can't run it with 4 ranks - I get

Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000010, n=4, ranks=0x555e6041cdb0, new_group=0x7ffefa6e5d20) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

2, 8, 16 ranks all work well with --backend=CPU ; 4 ranks also worked fine with 22.0.0

[stgeke]

Ok, looks like there was a bug in v22.0.
Can you please share a full backtrace of the 4 MPI rank case which fails.

[joneuhauser]

The log does not show more than this. Usually, if nekrs crashes, a backtrace is shown. How can I collect more information in this case?

relevant section of the log

 mesh metrics:
 GLL grid spacing min/max    : 5.00E-02 5.00E-01
 scaled Jacobian  min/max/avg: 1.00E+00 1.00E+00 1.00E+00
 aspect ratio     min/max/avg: 3.00E+00 1.00E+01 6.50E+00

================= INITIAL CONDITION ====================
U   min/max: 0 1  0 0  0 0
P   min/max: 0 0
S   min/max: 0 0
================= ELLIPTIC SETUP SCALAR00 ===============
bID 1 -> bcType codedFixedValue
bID 2 -> bcType zeroGradient
bID 3 -> bcType zeroGradient
timing gs: 7.15e-06s 1.22e-05s 1.55e-05s 1.57e-05s 1.69e-05s 
used config: pw (MPI min/max/avg: 2.86e-06s 3.10e-06s 3.04e-06s / avg bi-bw: 2.2GB/s/rank)
building Jacobi preconditioner ... done (6.48499e-05s)
done (0.00314403s)
================ ELLIPTIC SETUP VELOCITY ================
bID 1 -> bcType codedFixedValue
bID 2 -> bcType fixedGradient
bID 3 -> bcType zeroValue
timing gs: 1.03e-05s 2.38e-05s 2.31e-05s 2.34e-05s 2.41e-05s 
used config: pw (MPI min/max/avg: 3.10e-06s 4.05e-06s 3.64e-06s / avg bi-bw: 3.9GB/s/rank)
building Jacobi preconditioner ... done (6.10352e-05s)
done (0.00150752s)
================ ELLIPTIC SETUP PRESSURE ================
timing gs: 5.72e-06s 1.00e-05s 1.34e-05s 1.34e-05s 1.41e-05s 
used config: pw (MPI min/max/avg: 2.38e-06s 2.86e-06s 2.62e-06s / avg bi-bw: 1.7GB/s/rank)
building MG preconditioner ... 
============= BUILDING pMG2 ==================
timing gs: 5.72e-06s 9.78e-06s 1.26e-05s 1.26e-05s 1.34e-05s 
used config: pw (MPI min/max/avg: 2.15e-06s 2.38e-06s 2.26e-06s / avg bi-bw: 1.0GB/s/rank)
timing gs: 8.58e-06s 1.19e-05s 1.72e-05s 1.72e-05s 1.81e-05s 
used config: pw (MPI min/max/avg: 2.15e-06s 2.86e-06s 2.50e-06s / avg bi-bw: 1.5GB/s/rank)
estimating maxEigenvalue ... 1.35528 done (0.00135732s)
============= BUILDING COARSE pMG1 ==================
computing geometric factors ... J [0.00416666,0.0166667] done (6.55651e-05s)
meshParallelGatherScatterSetup N=1
done (6.7234e-05s)
timing gs: 4.29e-06s 7.87e-06s 9.78e-06s 9.78e-06s 1.05e-05s 
used config: pw (MPI min/max/avg: 1.91e-06s 2.15e-06s 2.09e-06s / avg bi-bw: 0.3GB/s/rank)
building FEM matrix ... done.
done (0.000266552s)
setup FEM solver ...Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000020, n=4, ranks=0x5581646291b0, new_group=0x7fff981dcca0) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000018, n=4, ranks=0x5611293724b0, new_group=0x7ffd25401aa0) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

[cli_1]: aborting job:
Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000018, n=4, ranks=0x5611293724b0, new_group=0x7ffd25401aa0) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000010, n=4, ranks=0x55aafe4220a0, new_group=0x7ffcecbacb60) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

[cli_2]: aborting job:
Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000010, n=4, ranks=0x55aafe4220a0, new_group=0x7ffcecbacb60) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000010, n=4, ranks=0x563819b22810, new_group=0x7ffcac1117c0) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

[cli_3]: aborting job:
Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000010, n=4, ranks=0x563819b22810, new_group=0x7ffcac1117c0) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

[cli_0]: aborting job:
Fatal error in PMPI_Group_incl:
Invalid rank, error stack:
PMPI_Group_incl(172).............: MPI_Group_incl(group=0xc8000020, n=4, ranks=0x5581646291b0, new_group=0x7fff981dcca0) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 2, has value 0 which is also the value at index 0

[stgeke]

How do you average along the homogenous directions? I can think of two ways: 1) planar_avg (works only for extruded or tensor-product meshes) or 2) pointInterpolation.

[joneuhauser]

Do you prefer a warning instead of an error?

I'd prefer a fix, but in general, warning about unmet assumptions in the code seems to be prudent.

Unfortunately there is no way to adjust ldimt automatically. We're working on a new native field file writer (instead of recycling the nek5000's implementation) to avoid this restriction.

That sounds promising! Is this the point where the fieldOffsets for velocity, temperature and scalars will actually be different? I'm wondering if I need to design my user code such that everything still works if this changes in the future.

Also, I think it would be useful to have an option (per scalar) that the other scalars, not only temperature, are also solved on the tmesh. One obvious application would be to simulate different thermal boundary conditions (where temperature is a passive scalar) at the same time in a CHT mesh. Is this something you have on the roadmap, or that I should look into?

[stgeke]

• At the moment the fieldOffset is same across all fields. However that might change in the future, when we support a different polynomialOrder per scalar.
• Supporting more than just one t-mesh scalar is of low-priority - at least for us.

[stgeke]

Should be fixed in latest next.

[joneuhauser]

Yes, the averaging procedure for the gabls case is very similarly to what I've written for my CHT case, thanks for the confirmation!

Just to confirm:

1. Do you prefer a warning instead of an error? Unfortunately there is no way to adjust ldimt automatically. We're working on a new native field file writer (instead of recycling the nek5000's implementation) to avoid this restriction.
2. Right let me fix this.

Those two are not fixed yet, right?

[stgeke]

I fixed (2) as well.