fluids: Add diffusive flux projection to Navier-Stokes

[jrwrigh]

Adds projection of the divergence of diffusive flux ($F^{diff}_{i,i}$) to the strong residual.
The projection is done by a lumped mass matrix multiplied against $\int_\Omega v_{,i} F^{diff}_i d \Omega$ . Note that there should be a boundary integral term here for consistencies ($\int_\Gamma v \cdot F^{diff}_i \hat n_i d\Gamma$).

TODOs:

• Add into Jacobian (currently just zeroed, not sure if it can/should be in the Jacobian) Shouldn't be done, see discussion below.
• Add Boundary Integral

[jedbrown]

Can you tighten the tolerance on the linear solve and test the reduction in nonlinear residual per Newton step? I think including this stabilization term shouldn't greatly change the nonlinear convergence, and thus I would leave it out, especially so long as we aren't using a very strong preconditioner. It can't be included in the assembled Jacobian without significant cost/loss of sparsity (if we were to assemble the exact Jacobian) and I think that definitely won't pay off. It may be good to document what approximations are made in the Jacobian. We don't differentiate through $\boldsymbol\tau$ and we don't include this term.

[jrwrigh]

Can you tighten the tolerance on the linear solve and test the reduction in nonlinear residual per Newton step? I think including this stabilization term shouldn't greatly change the nonlinear convergence, and thus I would leave it out, especially so long as we aren't using a very strong preconditioner.

Sure. Just do an A/B test of zeroing the diffusive flux and not?

It can't be included in the assembled Jacobian without significant cost/loss of sparsity (if we were to assemble the exact Jacobian) and I think that definitely won't pay off.

Yeah, that's the conclusion I was leaning towards.

It may be good to document what approximations are made in the Jacobian. We don't differentiate through and we don't include this term.

I can definitely do that, probably also include a blurb about the diffusive flux projection.

[jedbrown]

Yeah, like -ksp_rtol 1e-5 -snes_monitor -ksp_converged_reason with and without diffusive flux in the stabilization term. If we see a similar residual history, then we can be comfortable about not including it in the linear solve.

[jrwrigh]

I've added a flag -use_fluxproj_jac to turn the addition on and off. Interestingly, including the diffusive flux makes the convergence significantly worse. Using:

./navierstokes -options_file channel.yaml -{ts,snes}_monitor -{snes,ksp}_converged_reason -primitive -ts_max_steps 5 -stab supg -ksp_rtol 1e-5 -use_fluxproj_jac false

it converges within 4 non-linear iterations per timestep. But using -use_fluxproj_jac true causes it to take 50 before it reaches it's maximum non-linear iterations. I could bump it up more, but it hasn't even reduced the residual by an order of magintude (reduced from 2.07e-2 to 2.2e-3).

-use_fluxproj_jac false:

0 TS dt 5e-06 time 0.
    0 SNES Function norm 2.074482031710e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 7
    1 SNES Function norm 4.820634232613e-04 
    Linear solve converged due to CONVERGED_RTOL iterations 7
    2 SNES Function norm 1.369035192689e-05 
    Linear solve converged due to CONVERGED_RTOL iterations 6
    3 SNES Function norm 4.285282034630e-07 
  Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 3

use_fluxproj_jac true:

0 TS dt 5e-06 time 0.
    0 SNES Function norm 2.074482031710e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 22
    1 SNES Function norm 2.011876301805e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 22
    2 SNES Function norm 1.884363105444e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 23
    3 SNES Function norm 1.766406980818e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 23
    4 SNES Function norm 1.657118594118e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 23
    5 SNES Function norm 1.555844670806e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 23
    6 SNES Function norm 1.461973574009e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 23
    7 SNES Function norm 1.374940027875e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 22
    8 SNES Function norm 1.294223873158e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 20
    9 SNES Function norm 1.219349738033e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 22
   10 SNES Function norm 1.149874867108e-02 
    Linear solve converged due to CONVERGED_RTOL iterations 22
....

[jedbrown]

That sounds like a bug in -use_fluxproj_jac true. Can you compare with and without flux projection? I.e., leaving it out in the nonlinear residual versus adding it in. Keep the same Jacobian in both cases.

[jrwrigh]

Is it possible that it actually should be zero? Thinking from the integral form of the SUPG term:
$$\int_\overline{\Omega} \mathcal{P}(v) \tau (U_{,t} - S + F^{adv} + F^{diff}) d\overline{\Omega}$$.

We currently do take the derivative of $U_{,t}$, $S$, and $F^{adv}$, but don't for $\tau$ or for $\mathcal{P}(v)$. Thus the derivative of the term $\mathcal{P}(v) \tau F^{diff}$ in the integral is zero. Only thing I'm hazy on is whether it is true that we don't take the derivative of $\mathcal{P}(v)$.

[jedbrown]

Note that $F^{\text{adv}}$ and $F^{\text{diff}}$ are not fluxes (as the usual notation would imply), but the divergence of flux. I think we usually use $\mathcal L^{\text{adv}} Y$ and $\mathcal L^{\text{diff}} Y$ for those objects.

Yes, it's true that we "freeze" $L^{\text{adv}}$ as applied to the test function. I think this would not be especially hard to "fix" by computing this derivative exactly, but I wouldn't expect it to have a huge impact on nonlinear convergence.

[jrwrigh]

Note that $F^{\text{adv}}$ and $F^{\text{diff}}$ are not fluxes (as the usual notation would imply), but the divergence of flux. I think we usually use $\mathcal L^{\text{adv}} Y$ and $\mathcal L^{\text{diff}} Y$ for those objects.

Yep, you're right.

Yes, it's true that we "freeze" $L^{\text{adv}}$ as applied to the test function. I think this would not be especially hard to "fix" by computing this derivative exactly, but I wouldn't expect it to have a huge impact on nonlinear convergence.

I doubt that it'd have a significant impact either. Given that, would you agree that the behavior I saw with adding the divergence of diffusive flux is correct though? ie. we should leave it out of the jacobian (or more explicitly set it to zero in the function).

[jedbrown]

This branch had conflicts with main and will need a rebase or merge. Do you think it's close? We don't need this term in the Jacobian.

[jrwrigh]

The primary thing left here is to actually validate it against Ken's original paper. I had been trying to do that with the channel example, but ran into issues. There are two routes to doing that:

1. Periodic BCs and a non-equispaced grid
2. Inflow/Outflow BCs with an equispaced grid

Route 1 is stalled as I can't manipulate the grid with periodic BCs. Or at least I can't do that with the current methodology.

Route 2 is stalled in two ways. The channel inflow doesn't have a jacobian QF defined. I tried to run it with snes_fd_color, but that was unstable. I tried to quickly write up an inflow jacobian QF, but there was a bug in it such that it wouldn't converge. That work is local right now, but I can push it up if you want to take a look.

Running full production cases became priority, so I put this on hold until Route 1 becomes feasible, as getting the periodic BCs is necessary in the future anyways.

I also have some local work for getting the boundary integral term working. That could be another PR though if that's preferred.

[jedbrown]

I think periodic non-equally spaced should work, we just need to map the cell geometry and won't be able to evaluate surface integrals (for now). I'll try to fix the general periodic case in PETSc.

[jrwrigh]

I think periodic non-equally spaced should work, we just need to map the cell geometry and won't be able to evaluate surface integrals (for now).

Do you have a general idea of what needs to be done mapping the cell geometry?

[jedbrown]

Can you rebase (it's needed anyway) and leave a stub for what you've tried?

[jrwrigh]

Yep, working on that now.

[jrwrigh]

Branch is rebased now. It looks like I never commited my changes to get an non-equispaced grid. I can try to throw that together real quick.

[jrwrigh]

Actually I did commit them, but was looking at the wrong file history to find them. They're now pushed to jrwrigh/diff_flux_channelwork. That includes the buggy inflow jacobian QF and the mesh modification. (The inflow jacobian doesn't get used for the periodic BCs, so the changes are independent).

[jedbrown]

build errors? What should I run to test your attempt at mesh modification? (Let's ignore inflow/outflow for now.)

In file included from /home/jed/src/libCEED/examples/fluids/problems/channel.c:12:
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h: In function ‘Channel_Inflow’:
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:145:22: warning: implicit declaration of function ‘StateFromQi’; did you mean ‘StateFromY’? [-Wimplicit-function-declaration]
  145 |     State s_inside = StateFromQi(gas, q_inside, x);
      |                      ^~~~~~~~~~~
      |                      StateFromY
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:145:22: error: invalid initializer
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:125:20: warning: unused variable ‘gamma’ [-Wunused-variable]
  125 |   const CeedScalar gamma       = cp / cv;
      |                    ^~~~~
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h: In function ‘Channel_Inflow_Jacobian’:
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:203:16: warning: implicit declaration of function ‘StateFromQi_fwd’; did you mean ‘StateFromY_fwd’? [-Wimplicit-function-declaration]
  203 |     State ds = StateFromQi_fwd(gas, s, dqi, x_i, dx_i);
      |                ^~~~~~~~~~~~~~~
      |                StateFromY_fwd
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:203:16: error: invalid initializer
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:196:22: warning: unused variable ‘wdetJb’ [-Wunused-variable]
  196 |     const CeedScalar wdetJb  = (implicit ? -1. : 1.) * q_data_sur[0][i];
      |                      ^~~~~~
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:189:20: warning: unused variable ‘gamma’ [-Wunused-variable]
  189 |   const CeedScalar gamma             = cp / cv;
      |                    ^~~~~
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h: In function ‘Channel_Outflow’:
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:256:21: error: invalid initializer
  256 |     const State s = StateFromQi(&context->newtonian_ctx, qi, x_i);
      |                     ^~~~~~~~~~~
In file included from /home/jed/src/libCEED/include/ceed/ceed.h:72,
                 from /home/jed/src/libCEED/include/ceed.h:1,
                 from /home/jed/src/libCEED/examples/fluids/problems/../navierstokes.h:11,
                 from /home/jed/src/libCEED/examples/fluids/problems/channel.c:11:
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h: At top level:
/home/jed/src/libCEED/examples/fluids/problems/../qfunctions/channel.h:232:16: warning: ‘Channel_Outflow’ defined but not used [-Wunused-function]
  232 | CEED_QFUNCTION(Channel_Outflow)(void *ctx, CeedInt Q,
      |                ^~~~~~~~~~~~~~~
/home/jed/src/libCEED/include/ceed/types.h:48:34: note: in definition of macro ‘CEED_QFUNCTION’
   48 |   CEED_QFUNCTION_ATTR static int name
      |                                  ^~~~

[jrwrigh]

Sorry about that. Should build fine now. I've moved the inflow/outflow stuff into a different (local) branch for right now.

[jrwrigh]

I'm having trouble getting a convergent channel simulation. My previous copy of the branch is running fine, so I'm guessing I missed something in the rebasing that didn't get caught by the tests. I'll let you know when the branch is actually ready to look at.

[jrwrigh]

Turned out to just be a -primitive vs -state_var primitive issue. Rebase was fine. This is the code that I have working on jrwrigh/diff_flux_channelwork

./navierstokes -options_file channel.yaml -{ts,snes}_monitor -{snes,ksp}_converged_reason -state_var primitive -dm_plex_box_faces 20,20,1 -stab supg

[jedbrown]

This makes it look like the boundary integral might have been applied to the bottom face, but not top face?

Or is that only an effect of the grid mapping?

I think we also need to DMGetCellCoordinatesLocal and apply the same mapping. I'll tinker with the code if that doesn't solve it for you.

[jrwrigh]

There isn't an boundary integral applied anywhere (otherwise getting the element restriction would fail). It's a result of the mesh mapping:

Left is with the mapping, right is without.

I'll try with the DMGetCellCoordinatesLocal and see what that does.

[jrwrigh]

Yep, that did it:

Note that the values returned by DMGetBoundingBox doesn't include the last layer of elements on the periodic face. I understand that the faces of those elements technically "don't exist" since they're technically the same as their periodic counter part, but I feel like the rest of the element still counts. If we're thinking about it from the perspective of the domain mapping onto a torus, I want the full circumference of the torus, not the circumference of the torus minus some arc-length.

Anywho, I'll try and get validation results soon.

[knepley]

@jrwrigh You are correct that the bounding box calculation was wrong for periodic things. I believe it is now fixed.

[jrwrigh]

@jrwrigh You are correct that the bounding box calculation was wrong for periodic things. I believe it is now fixed.

Ok, cool. I'm currently on main, 6fcee9a (from August 9th). Do you know of a specific MR/commit that should've fixed it? Quick search didn't show anything, but I might not have the right keywords to make it pop up. I'll pull down the latest PETSc and rebuild it tomorrow to make sure I haven't missed anything.

[knepley]

Dang, you are right. I fixed it in this unfortunate branch that is taking forever: https://gitlab.com/danofinn/petsc/-/commits/danfinn/feature-swarm-poisson-test

[jedbrown]

@jrwrigh I don't think it's necessary to actually do the surface integrals, just to zero the wall nodes. The reason comes from the identity that $\int_{\Omega} \nabla\cdot F = \int_{\partial\Omega} F \cdot n$, which holds also at the element scale. In solid mechanics, these are called "reaction forces", and are a more accurate way to evaluate forces than to literally do the surface integral.