Planar averaging for unstructured mesh

[moviegroovie]

Hello,

I am interested in performing planar averaging of fields, similar to what is done by the planarAvg function in the gabls example. However, I'd like to be able to perform planar averaging on custom meshes that may not be structured. For example, for a flow-over-cylinder geometry, my mesh would be unstructured in the cross-sectional plane, and I would like to take averages in the direction of the cylinder axis. Firstly, is my understanding correct that the planarAvg function assumes that the mesh is made with genbox and a .box file? Secondly, is there already a way to perform planar averaging in NekRS for custom/unstructured meshes? If not, is there a simple way to modify the existing planar averaging function to obtain this desired functionality?

Thank you,
Hoyean

[yslan]

The built in planar average relies on the topology of the mesh and it doesn't work for unstructured grid.

Two ways for you.

1. Run a dummy NekRS simulation with a structured cylinder mesh so it's a 2D mesh then extruded into the main axis. Then, load your data via gfldr or via restartFrom = "a.fld + int" (next branch). Since the dummy simulation uses structural grid, you can call planar average.
2. Use findpt to do interpolation your self. Define your own grid for the numerical integral, findpt interpolate the solution onto those grid, then you can have full control.
   

   nekRS/examples/gabls1/gabls.udf

   Lines 212 to 243 in 55d6414

   	static pointInterpolation_t *interpolator = nullptr;
   	if (!interpolator) {
   	const auto np = (platform->comm.mpiRank() == 0) ? NUMBER_ELEMENTS_Y * mesh->N : 0;
   	xp.resize(np);
   	yp.resize(np);
   	zp.resize(np);
   	const auto yMin = platform->linAlg->min(mesh->Nlocal, mesh->o_y, platform->comm.mpiComm());
   	const auto yMax = platform->linAlg->max(mesh->Nlocal, mesh->o_y, platform->comm.mpiComm());
   	if (np) {
   	const auto x0 = 0.0;
   	const auto z0 = 0.0;
   	for (int i = 0; i < np; i++) {
   	xp[i] = x0;
   	yp[i] = static_cast<dfloat>(i)/(np-1) * (yMax-yMin);
   	zp[i] = z0;
   	}
   	}
   	interpolator = new pointInterpolation_t(mesh, platform->comm.mpiComm());
   	interpolator->setPoints(xp, yp, zp);
   	interpolator->find(pointInterpolation_t::VerbosityLevel::Detailed);
   	if (platform->comm.mpiRank() == 0) {
   	std::ofstream f(fname, std::ios::app);
   	f << R"(# time y <u> <w> <T> <dudy> <dwdy> <dTdy> <u'T'> <w'T'> <w'w'T'>)" << std::endl;
   	f.close();
   	}
   	}

We are preparing a more user friendly plugin as an interface for 2. so everyone can extract the data on their favored way easier.

[moviegroovie]

Got it, thank you!

[moviegroovie]

Hello,

Can the planarAvg routines be modified for mesh topology being read from gmsh? The topology data exists (if you mean connnectivity by topology). We would be happy to give it a shot if you provide explanation on how the topology is used. Thank you very much for your insights.

[yslan]

Nek5000 uses "lexicographical ordering" for a single box or an extruded mesh.

• Box mesh
  Say we have a box mesh with nelx, nely and nelz elements in x y and z direction.
  The elements are ordered in this way:

  integer ex, ey, ez, e
  
  ! go through all elements 
  do ez = 1, nelz
  do ey = 1, nely
  do ex = 1, nelx
    eg = (ez-1) * nelx * nely + (ey-1) * nelx + ex  ! global element ordering 
  enddo
  enddo
  enddo
  
  ! go through locally stored elements within this MPI rank
  do e=1,nelt
    eg = lglel(e)
    call get_exyz(ex,ey,ez,eg,nelx,nely,nelz) ! from eg to (ex,ey,ez)
  enddo
  

• Extruded mesh
  Suppose we have a 2D mesh on x-y plane with nelxy elements. The 3D mesh is extruded in z direction using nelz elements. Then the elements are ordered

  integer exy, ez, e, etmp
  
  ! go through all elements 
  do ez = 1, nelz
  do exy = 1, nelxy
    eg = (ez-1) * nelxy + exy  ! global element ordering 
  enddo
  enddo
  
  ! go through locally stored elements within this MPI rank
  do e=1,nelt 
    eg = lglel(e)
    call get_exyz(exy,etmp,ez,eg,nelxy,1,nelz) ! from eg to (exy,ez)
  enddo
  

Same for NekRS except the index are now 0-base rather than 1-base.

planarAvg uses this known ordering to do line integrals.
The key is sum. Because we know how the grid points are ordered, we can assign integer to grid points such that all points along the same "fiber" share the same integer. Then, we can do gather-scatter, or thinking of Allreduce under that integer tag.