Clarification on opSEM usage

[foreverkon]

Dear all,

I am relatively new to NekRS and have been studying the examples in the next branch. I encountered some confusion regarding the usage of functions within opSEM.

Here are my understanding:

1. Functions prefixed with strong (e.g., strongCurl, strongGrad) compute the strong form of the derivative. By default, the avg parameter is true, so they return C0 continuous fields (physical values).
2. Functions without the prefix compute the weak form. These return discontinuous, element-wise results that are weighted by the $Jw$.
3. The runAvg function first calls startFinish (gather-scatter option) to sum contributions at shared nodes, and then multiplies by o_invAJw (the inverse of the assembled global mass matrix) to recover the consistent strong form values.

But I am still confused about the following code snippets.
In the examples/tgv, the enstrophy is computed as follows:

const auto enstrophy = [&] () {
    const dfloat scale = 0.5 / mesh->volume;
    auto o_curlU = opSEM::strongCurl(mesh, nrs->fieldOffset, nrs->fluid->o_U);

    static poolDeviceMemory<dfloat> o_invJw;
    if (!o_invJw.isInitialized()) {
        o_invJw.resize(mesh->Nlocal);
        platform->linAlg->adyz(mesh->Nlocal, 1.0, mesh->o_Jw, o_invJw);
    }

    platform->linAlg->axmyMany(mesh->Nlocal, mesh->dim, nrs->fieldOffset, 0, 1.0, o_invJw, o_curlU);

    poolDeviceMemory<dfloat> o_magCurlU(nrs->fieldOffset);
    platform->linAlg->magSqrVector(mesh->Nlocal, nrs->fieldOffset, o_curlU, o_magCurlU);

    return scale * platform->linAlg->innerProd(mesh->Nlocal,
                        o_magCurlU,
                        mesh->o_Jw,
                        platform->comm.mpiComm(),
                        0);
    }();

I am confused about the usage of o_invJw. If strongCurl (with default avg=true) already returns the physical vorticity field $\boldsymbol{\omega}$, why is it necessary to multiply by o_invJw?
Mathematically, this seems to compute $\sum (\boldsymbol{\omega}/Jw)^2 Jw$, whereas the definition of enstrophy is
$\sum \boldsymbol{\omega}^2 Jw$. This logic would only seem consistent if strongCurl returned the weighted weak form (i.e., if avg was false).

I noticed a somewhat related pattern in the HIT example for dissipation calculation:

const auto epsVolAvg = [&]() {
    dfloat mue, rho;
    platform->options.getArgs("FLUID VISCOSITY", mue);
    platform->options.getArgs("FLUID DENSITY", rho);
    const auto nu = mue / rho;

    auto o_gradU = opSEM::strongGradVec(mesh, nrs->fieldOffset, nrs->fluid->o_U);

    static poolDeviceMemory<dfloat> o_invAJwJw;
    if (!o_invAJwJw.isInitialized()) {
        o_invAJwJw.resize(mesh->Nlocal);
        platform->linAlg->axmyz(mesh->Nlocal, 1.0, mesh->o_invAJw, mesh->o_Jw, o_invAJwJw);
    }
    platform->linAlg->axmyMany(mesh->Nlocal, mesh->dim * mesh->dim, nrs->fieldOffset, 0, 1.0, mesh->o_invAJw, o_gradU);
    nrs->qqt->startFinish("+", o_gradU, nrs->fieldOffset);

    poolDeviceMemory<dfloat> o_eps(mesh->Nlocal);
    platform->linAlg->magSqrTensor(mesh->Nlocal, nrs->fieldOffset, o_gradU, o_eps);
    platform->linAlg->scale(mesh->Nlocal, nu, o_eps);
    return avg(o_eps);
    }();

Here, ignoring o_invAJwJw (which seems unused), after strongGradVec, the code multiplies by o_invAJw and explicitly calls startFinish. Since strongGradVec presumably already performs gather-scatter, this looks like an additional smoothing/filtering step. I am trying to understand the motivation behind this.

My questions are:

1. Is my understanding of the strong vs. non-strong functions correct?
2. In the examples above, why are these additional scalings (o_invJw in TGV, o_invAJw in HIT) applied after the strong derivative calculation? What is the key distinction between using o_invJw vs o_invAJw in these post-processing contexts?
3. If my goal is simply to compute the time-average of the velocity gradient, is it correct to pass the result of strongGradVec directly to the tavg module without extra scaling?

I really appreciate any clarifications you can provide on these points. Thank you very much!

[stgeke]

opSEM::strongCurl computes the curl using strong derivatives. However, the resulting field can still weighted by Jw (in case it's not smoother). By default, it's averaged across element boundaries to ensure C⁰ continuity and not weighted by Jw.

[yslan]

Please see this documentation for opSEM. Let me know if you need more clarification.
https://nekrs.readthedocs.io/en/latest/other_resources.html#sem-operators-opsem

Regarding o_Jw and o_invAJw under mesh_t. We have a table here:
https://nekrs.readthedocs.io/en/latest/other_resources.html#mesh-data

"Jw" stands for Jacobian times quadrature weight, which is the entries of the diagonal lumped mass matrix
"invAJw" is the inverse of the asssembled lumped mass matrix.
They are not the inverse to each other. If you multiple them, it yields the inverse of the multiplicity

[foreverkon]

Thanks for the explanation, that really helps.
I hadn't noticed the recent updates to the documentation until now.

But I’m still a bit confused about two specific points:

1. Why do we need o_invJw and o_invAJw in the enstrophy and dissipation examples mentioned above?
2. Regarding the docs (https://nekrs.readthedocs.io/en/latest/other_resources.html#sem-operators-opsem), it says strongGrad (when avg=true) calculates $[ G D_{x} u , G D_{y} u , G D_{z} u ]$, but I feel like it should be $[ GBD_{x} u , GBD_{y} u , GBD_{z} u ]$. Could you confirm?

[yslan]

1. There is a missing "false" in tgv example. This has been fixed in the private development branch but hasn't been put to next.

    auto o_curlU = opSEM::strongCurl(mesh, nrs->fieldOffset, nrs->fluid->o_U, false);
   

   It's simply to demonstrate how to recover a unweighted term from a weighted strong form. It's not $C^0$ but it's not that important for enstrophy since there is a summation later.

2. strongGrad with avg=true is $GDu$ without weights. See the relevant code:
   kernel:
   

   nekRS/src/core/kernels/gradientVolumeHex3D.okl

   Lines 91 to 94 in b5b876d

   	const dlong id = e * p_Np + k * p_Nq * p_Nq + j * p_Nq + i;
   	gradP[id + 0 * offset] = JW * (drdx * dpdr + dsdx * dpds + dtdx * dpdt);
   	gradP[id + 1 * offset] = JW * (drdy * dpdr + dsdy * dpds + dtdy * dpdt);
   	gradP[id + 2 * offset] = JW * (drdz * dpdr + dsdz * dpds + dtdz * dpdt);

   
   opSEM::strongGrad
   

   nekRS/src/core/opSEM.cpp

   Lines 34 to 44 in b5b876d

   	void strongGrad(mesh_t *mesh, dlong offset, const occa::memory &o_in, occa::memory &o_out, bool avg)
   	{
   	launchKernel(section + "gradientVolume" + suffix,
   	mesh->Nelements,
   	mesh->o_vgeo,
   	mesh->o_D,
   	offset,
   	o_in,
   	o_out);
   	if (avg) runAvg(mesh, o_out, offset, mesh->dim);
   	}

   
   avg=true, startFinish does the gather-scatter. o_invAJw contains $B^{-1}$ and $M^{-1}$.
   

   nekRS/src/core/opSEM.cpp

   Lines 7 to 11 in b5b876d

   	static void runAvg(mesh_t *mesh, occa::memory& o_U, dlong offset, int nFields)
   	{
   	oogs::startFinish(o_U, nFields, offset, ogsDfloat, ogsAdd, mesh->oogs);
   	platform->linAlg->axmyMany(mesh->Nlocal, nFields, offset, 0, 1.0, mesh->o_invAJw, o_U);
   	}

[foreverkon]

Understood. Thanks.