Anisotropic implicit forcing

[makrandak]

Hello,
I am on the latest commit of next branch. I want to say first of all the feature of userVelocityImplicitLinearTerm is very helpful as it makes sure that adding forcing does not ruin the accuracy of linear multistep method.

I am a bit confused by the setup of nrs->userVelocityImplicitLinearTerm though. Is my understanding correct that the code here assumes that the forcing will be the same for all three components of velocity.

Assuming we need the flexibility of having different forcing coefficients in different velocity components, in tombo.cpp, shouldn't lamba1 allocate memory for all three components like
auto lambda1 = platform->o_memPool.reserve<dfloat>(nrs->NVfields * nrs->fieldOffset); instead of auto lambda1 = platform->o_memPool.reserve<dfloat>(mesh->Nlocal);,
then if the user sets different coefficients for different components of velocity it can be read like this,

auto o_implicitLT = nrs->userVelocityImplicitLinearTerm(time);
const auto o_implicitLT_x = o_implicitLT.slice(0 * nrs->fieldOffset);
const auto o_implicitLT_y = o_implicitLT.slice(1 * nrs->fieldOffset);
const auto o_implicitLT_z = o_implicitLT.slice(2 * nrs->fieldOffset);
platform->linAlg
          ->axpbyz(mesh->Nlocal, nrs->g0 / nrs->dt[0], nrs->o_rho, 1.0, o_implicitLT, o_lambda1_x);
platform->linAlg
          ->axpbyz(mesh->Nlocal, nrs->g0 / nrs->dt[0], nrs->o_rho, 1.0, o_implicitLT, o_lambda1_y);
platform->linAlg
          ->axpbyz(mesh->Nlocal, nrs->g0 / nrs->dt[0], nrs->o_rho, 1.0, o_implicitLT, o_lambda1_z);

and subsequently used in the solver like this,

nrs->uSolver->solve(o_lambda0, o_lambda1_x, o_rhsX, o_U.slice(0 * nrs->fieldOffset));
nrs->vSolver->solve(o_lambda0, o_lambda1_y, o_rhsY, o_U.slice(1 * nrs->fieldOffset));
nrs->wSolver->solve(o_lambda0, o_lambda1_z, o_rhsZ, o_U.slice(2 * nrs->fieldOffset));

And for example in such a case the user could then write an implicit forcing like this in the .udf file.

occa::memory implicitForcing(dfloat time)
{
  auto mesh = nrs->meshV;
  std::vector<dfloat> AcompRead(3);
  platform->par->extract("casedata", "A11", AcompRead[0]);
  platform->par->extract("casedata", "A22", AcompRead[1]);
  platform->par->extract("casedata", "A33", AcompRead[2]);


  poolDeviceMemory<dfloat> o_F(nrs->NVfields * nrs->fieldOffset);
  
  for(int i = 0; i < mesh->dim; i++) {
    auto o_inlt = o_F.slice(i * nrs->fieldOffset, mesh->Nlocal);
    platform->linAlg->fill(mesh->Nlocal, AcompRead[i], o_inlt);
  }

  return o_F;
}

with
nrs->userVelocityImplicitLinearTerm = &implicitForcing; in UDF_Setup().

Let me know if this can already be done in the current version of the code or if I am misinterpreting the code.

I am very appreciative of your help and insight. Also, if this is useful, let me know if I can code it and contribute to the repo. Looking forward to hearing from you.

With sincere regards,
Makrand

[stgeke]

The implicit treatment will potentially improve stability but not accuracy. Currently, we do not support having different coefficients for each velocity component. Could you please specify the application you have in mind?

[makrandak]

The implicit treatment will potentially improve stability but not accuracy.

If the explicit source term is computed just at the last time state (t^{n}), which is my understanding if treating explicitly in the code, accuracy is lost. Maybe I am misunderstanding on whether explicit source terms are extrapolated to {t^n+1} state or not (for BDF2, two previous time states would be used to extrapolate to t^{n+1}). For BDF schemes an approximation of the state at the current time (t^{n+1}) preserves the BDF accuracy, which the implicit treatment provides.

Currently, we do not support having different coefficients for each velocity component. Could you please specify the application you have in mind?

Thank you, I work on probing turbulence (using DNS and LES) using various types of forcing in various scenarios to improve RANS models. It can also be used for synthetic turbulence, see for e.g. https://pubs.aip.org/aip/pof/article/27/3/035115/259740/Anisotropic-linear-forcing-for-synthetic
There are more examples where this is useful, which we can maybe discuss over email.

If you comment on whether the coding strategy I suggested above to accommodate anisotropic forcing, and if I have understood the code correctly. I can incorporate the changes in my fork. Although, it would be great if it is included in some future commit.

Looking forward to your thoughts.

With sincere regards,
Makrand

[stgeke]

We're using a high-order extrapolation scheme that preserves the overall accuracy of the BDF scheme. However, the explicitly treated terms can limit stability. If the term is linear, we can avoid this issue by treating it implicitly. What you're suggesting will be sufficient, but note that the Helmholtz operator will also need to be changed.

[makrandak]

We're using a high-order extrapolation scheme that preserves the overall accuracy of the BDF scheme. However, the explicitly treated terms can limit stability.

Okay, then I agree that the accuracy would be the same.

If the term is linear, we can avoid this issue by treating it implicitly. What you're suggesting will be sufficient, but note that the Helmholtz operator will also need to be changed.

The terms are linear think of something like A_ij (u_j - <u_j>). Linear but different in three directions.
Can you please point me to where else the modifications I suggested also need to be changed in the code? I thought just having different o_lambda1 in x, y, z would work. If you point me to the place in code. I will study it.

Thanks for your help.

With sincere regards,
Makrand

[fischer1]

The problem with the addition of an arbitrary combination of linear terms is that you cannot guarantee that the block 3x3 system will be SPD and amenable to conjugate gradient iteration. For example, Coriolis forces are linear but they are completely skew symmetric. Combining them with the linear diffusion operator would lead to complex eigenvalues in that solver, rather than positive, which is what we have for the 3x3 block system now. The other thing to keep in mind is that the eigenvalues of the diffusion operator are large—they scale roughly as dt*nu/(dx_min)^2, which is one of the reasons to treat them implicitly. External forcing is generally determined by physical considerations and tend not to increase as the mesh is refined, so the need to be implicit is not as great. Nonetheless, stabilizing where possible is a worthy objective.

…

________________________________ From: Makrand Khanwale ***@***.***> Sent: Saturday, June 22, 2024 3:04 PM To: Nek5000/nekRS ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [Nek5000/nekRS] Anisotropic implicit forcing (Discussion #569) We're using a high-order extrapolation scheme that preserves the overall accuracy of the BDF scheme. However, the explicitly treated terms can limit stability. Okay, then I agree that the accuracy would be the same. If the term is linear, we can avoid this issue by treating it implicitly. What you're suggesting will be sufficient, but note that the Helmholtz operator will also need to be changed. The terms are linear think of something like A_ij (u_j - <u_j>). Linear but different in three directions. Can you please point me to where else the modifications I suggested also need to be changed in the code? I thought just having different o_lambda1 in x, y, z would work. If you point me to the place in code. I will study it. Thanks for your help. With sincere regards, Makrand — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/Nek5000/nekRS/discussions/569*discussioncomment-9848748__;Iw!!DZ3fjg!--Oi373xR6JwAcI0m-do_kipj4E97SB-CLzsuvgD3KshA0t5B1e67k-GZNAc2wYC7Uvcvj_hHH4TY_-sP5efpIkCqwQ$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACV35EN2O6GEPTLR5H3YDXDZIXKGVAVCNFSM6AAAAABJW7BNB6VHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4TQNBYG42DQ__;!!DZ3fjg!--Oi373xR6JwAcI0m-do_kipj4E97SB-CLzsuvgD3KshA0t5B1e67k-GZNAc2wYC7Uvcvj_hHH4TY_-sP5efw5GUbIo$>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[makrandak]

Hello Paul,

In general, I agree with you. However, we have all A_ij as positive coefficients (but anisotropic) in this case, with the forcing A_ij (u_j - <u_j>) so it would not destroy the conditioning of the resultant operator, and this support is not a very far-fetched generalization of the current implementation. Such forcing techniques allow us to probe how turbulence develops when forced anisotropically.

I am also flexible to just implement the extension in my fork, but I also wanted to check if my implementation strategy was missing something, which is one of the main objectives of this thread. If there are places where code needs modification and I am missing them (based on the implementation strategy I detailed above), that would be immensely valuable information.

I also thought it would be a valuable generalization as long as there is a warning to the user that the coefficients have to be positive.

Thank you for looking into this and for your help and time.

With sincere regards,
Makrand

[fischer1]

Hi Makrand, We should be able to support that. All you actually need is for the 3x3 A_ij matrices to be SPD, and there is an easy check for that (if we want to bother with that). When you refer to anisotropy, I'm assuming you mean in the linear, non-differentiated, coefficients, rather than in the coefficients that are in the stress tensor? Anisotropy in the stress tensor can also be supported, but that is a different question. hth, Paul

…

________________________________ From: Makrand Khanwale ***@***.***> Sent: Wednesday, June 26, 2024 6:21 PM To: Nek5000/nekRS ***@***.***> Cc: Fischer, Paul ***@***.***>; Comment ***@***.***> Subject: Re: [Nek5000/nekRS] Anisotropic implicit forcing (Discussion #569) Hello Paul, In general, I agree with you. However, we have all A_ij as positive coefficients (but anisotropic) in this case, with the forcing A_ij (u_j - <u_j>) so it would not destroy the conditioning of the resultant operator, and this support is not a very far-fetched generalization of the current implementation. Such forcing techniques allow us to probe how turbulence develops when forced anisotropically. I am also flexible to just implement the extension in my fork, but I also wanted to check if my implementation strategy was missing something, which is one of the main objectives of this thread. If there are places where code needs modification and I am missing them (based on the implementation strategy I detailed above), that would be immensely valuable information. I also thought it would be a valuable generalization as long as there is a warning to the user that the coefficients have to be positive. Thank you for looking into this and for your help and time. With sincere regards, Makrand — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/Nek5000/nekRS/discussions/569*discussioncomment-9887973__;Iw!!DZ3fjg!7qf7-t5q2n3eAJvF7llBLsJzYBIqC4bDRkvOUWtQO4DcfN3XNr-ieqLHnNTszs4dvwi09CCtBSB4sjASZ4er_purQjY$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACV35EJNJZ2T2CAABGHNPNTZJNEHLAVCNFSM6AAAAABJW7BNB6VHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4TQOBXHE3TG__;!!DZ3fjg!7qf7-t5q2n3eAJvF7llBLsJzYBIqC4bDRkvOUWtQO4DcfN3XNr-ieqLHnNTszs4dvwi09CCtBSB4sjASZ4er-VSDNak$>. You are receiving this because you commented.Message ID: ***@***.***>

[makrandak]

When you refer to anisotropy, I'm assuming you mean in the linear, non-differentiated, coefficients, rather than in the coefficients that are in the stress tensor?

Yes, these are not coefficients in the stress tensor, these are coefficients of the linear forcing (different in each direction), all positive.

Thanks for your help.

With sincere regards,
Makrand