Fluids - Compute component-wise error for NS solver

[LeilaGhaffari]

This is part of a course project in collaboration with @jnclement and @amandashack . It adds support for computing the error for each component and converts the errors to the other set of variables (primitive <-> conservative). It currently breaks for the other solvers though (advection and euler).

[LeilaGhaffari]

I think this branch was working for our purposes in the course project but is broken for the Advection2D and Euler solvers. A quick fix would be to define a convert_state QFunction for each of those solvers. We could, instead, refactor the Euler problems to use newtonian.c (could we do the same thing with Advection?). However, the branch would diverge from its original purpose.

[jedbrown]

Hmm, I think State{,Primitive,Conservative} need not be "Newtonian" things (they have nothing to do with viscosity). Their conversions assume ideal gasses. So we could promote that concept to something that Euler would also use. That's scope creep and I think it would also be fine to just disable your feature for Advection and Euler.

[LeilaGhaffari]

Oh, true! I was talking about Newtonian since that is the only solver that uses NewtonianIdealGasContext, IIRC. Then, I would prefer disabling the feature for Advection and Euler for now as we have decided to refactor Euler (to use Newtonian) anyway.

[jedbrown]

Sounds good to me. Newtonian is what matters for us (and at some point in the next few months, I think we'll add RANS).

[jedbrown]

Looks good. Perhaps a topic for another day, but I wonder if the velocity error norm should be defined using the pointwise relative error $\frac{\lVert \mathbf u_h - \mathbf u_* \rVert}{\lVert \mathbf u_* \rVert}$ (or absolute vector error) so we don't have anomalies for exactly or approximately axis-aligned flows. I think this (from Blasius) would show less than 1% error across all fields if measured in that way.

  Primitive variables:
        Component 0: 6.06303e-05 (24.4948, 404002.)
        Component 1: 0.00177814 (0.276762, 155.647)
        Component 2: 5.3239 (0.256163, 0.0481156)
        Component 3: -nan. (0., 0.)
        Component 4: 0.00184015 (2.93409, 1594.48)

We could do that in postprocessing by comparing absolute error as a vector to the true solution norm as a vector. (What's happening here is that velocity is very accurate in total, but not in exactly the same direction as the exact solution, thus we get a small y component error that is "big" compared to the exact y component that is almost perfectly aligned with the axis.)

[LeilaGhaffari]

@jedbrown , does this look good to you now?

[LeilaGhaffari]

@jedbrown , can we merge this branch if there are no other concerns?

[jedbrown]

Looks pretty close. Could you test or at least show sample output, perhaps using the Euler vortex?

[LeilaGhaffari]

Here is the output with the traveling vortex problem:

$ build/fluids-navierstokes -problem euler_vortex -degree 2 -q_extra 2 -dm_plex_box_faces 5,5,1 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 125,125,250 -dm_plex_dim 3 -mean_velocity 1.4,-2.,0 -bc_inflow 4,6 -bc_outflow 3,5 -bc_slip_z 1,2 -ts_dt 1e-7 -ts_rk_type 5bs -ts_rtol 1e-10 -ts_atol 1e-10 -ts_max_time 1

-- Navier-Stokes solver - libCEED + PETSc --
  MPI:
    Host Name                          : leila-ThinkPad-P53s
    Total ranks                        : 1
  Problem:
    Problem Name                       : euler_vortex
    Test Case                          : isentropic_vortex
    Background Velocity                : 1.400000,-2.000000,0.000000
    Stabilization                      : none
  libCEED:
    libCEED Backend                    : /cpu/self/opt/blocked
    libCEED Backend MemType            : host
  PETSc:
    Box Faces                          : 5,5,1
    DM MatType                         : aij
    DM VecType                         : standard
    Time Stepping Scheme               : explicit
  Mesh:
    Number of 1D Basis Nodes (P)       : 3
    Number of 1D Quadrature Points (Q) : 5
    Global DoFs                        : 1573
    Owned DoFs                         : 1573
    DoFs per node                      : 5
    Global nodes (DoFs / 5)            : 314
    Local nodes                        : 363
Time taken for solution (sec): 1.41848

Relative Error (absolute error norm, true solution norm):
  Conservative variables:
        Component 0: 0.0471095 (0.047109, 0.99999)
        Component 1: 0.364438 (0.505128, 1.38605)
        Component 2: 0.133585 (0.299979, 2.24561)
        Component 3: inf. (5.52411e-16, 0.)
        Component 4: 0.0467154 (0.264258, 5.65677)

Time integrator took 245 time steps to reach final time 1.00194

[jedbrown]

Why does it always say Conservative variables (regardless of -state_var primitive) and wouldn't we expect this to converge under refinement? Also, the domain shape (dm_plex_box_upper 125,125,250) relative to resolution dm_plex_box_faces 5,5,1 looks weird.

[LeilaGhaffari]

Why does it always say Conservative variables (regardless of -state_var primitive) ?

It's because traveling vortex doesn't support primitive variables yet.

[LeilaGhaffari]

wouldn't we expect this to converge under refinement?

I remember we didn't see convergence with -degree 2 in the summer. I only presented degree 1 in my talk.

[LeilaGhaffari]

Also, the domain shape (dm_plex_box_upper 125,125,250) relative to resolution dm_plex_box_faces 5,5,1 looks weird.

I took the arguments from the regression tests in navierstokes.c (and refined it) to show you the output. Oops, it is indeed weird!

[LeilaGhaffari]

Can we have one test of this functionality?

Do you mean a convergence plot?

[jedbrown]

Well, something that runs this code and checks that the errors are of expected size. (What if a bug caused the reported error to be zero always or the error became huge?) Ideally we would make a convergence plot for the isentropic vortex or channel flow, but I can handle that not being ready.

[jrwrigh]

Thanks for getting this done!

[jedbrown]

Cool. Squash or squash-merge?

[LeilaGhaffari]

Thanks, @jrwrigh for bearing with me. A lot has changed since the last time I contributed to this mini-appand. There are still a couple of remaining tasks before merging:

• Adding a test that Jed mentioned above
• Fixing the broken code for traveling vortex and 2d advection
• Fixing the bug in reporting the error (the computed error depends on the mpi rank)

[LeilaGhaffari]

Cool. Squash or squash-merge?

squash-merge is fine after I am done with the remaining tasks.

[LeilaGhaffari]

I am seeing a weird issue. The errors are not correct with rank 1.

$ build/fluids-navierstokes -problem euler_vortex -degree 3 -dm_plex_box_faces 1,1,2 -dm_plex_box_lower 0,0,0 -dm_plex
_box_upper 125,125,250 -dm_plex_dim 3 -units_meter 1e-4 -units_second 1e-4 -mean_velocity 1.4,-2.,0 -bc_inflow 4,6 -bc_outflow 3,5 -bc_slip_z 1,2 -vortex_strength 2 -ksp_atol 1e-4 -ksp_rtol 
1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 -implicit -dm_mat_preallocate_skip 0 -ts_type alpha


L2 Error:
  Conservative variables-Component 0: 0.000741784
  Conservative variables-Component 1: 0.000961908
  Conservative variables-Component 2: 0.00153698
  Conservative variables-Component 3: 0.
  Conservative variables-Component 4: 0.00399707
Time integrator CONVERGED_TIME on time step 50 with final time 0.05

$ mpiexec -n 4  build/fluids-navierstokes -problem euler_vortex -degree 3 -dm_plex_box_faces 1,1,2 -dm_plex_box_lo
wer 0,0,0 -dm_plex_box_upper 125,125,250 -dm_plex_dim 3 -units_meter 1e-4 -units_second 1e-4 -mean_velocity 1.4,-2.,0 -bc_inflow 4,6 -bc_outflow 3,5 -bc_slip_z 1,2 -vortex_strength 2 -ksp_at
ol 1e-4 -ksp_rtol 1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 -implicit -dm_mat_preallocate_skip 0 -ts_type alpha

L2 Error:
  Conservative variables-Component 0: 2.66349e-06
  Conservative variables-Component 1: 8.8622e-05
  Conservative variables-Component 2: 6.99578e-05
  Conservative variables-Component 3: 3.72699e-155
  Conservative variables-Component 4: 2.05213e-05
Time integrator CONVERGED_TIME on time step 50 with final time 0.05

$ mpiexec -n 8  build/fluids-navierstokes -problem euler_vortex -degree 3 -dm_plex_box_faces 1,1,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 125,125,250 -dm_plex_dim 3 -units_meter 1e-4 -units_second 1e-4 -mean_velocity 1.4,-2.,0 -bc_inflow 4,6 -bc_outflow 3,5 -bc_slip_z 1,2 -vortex_strength 2 -ksp_atol 1e-4 -ksp_rtol 1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 -implicit -dm_mat_preallocate_skip 0 -ts_type alpha

L2 Error:
  Conservative variables-Component 0: 2.66349e-06
  Conservative variables-Component 1: 8.8622e-05
  Conservative variables-Component 2: 6.99578e-05
  Conservative variables-Component 3: 4.56152e-155
  Conservative variables-Component 4: 2.05213e-05
Time integrator CONVERGED_TIME on time step 50 with final time 0.05

(same issue with advection2d and blasius)
It could be the logic in FixMultiplicity_NS() but I can't point at anything specific.

[LeilaGhaffari]

Testing with different arguments/problems, I see that the results could be different for n>1 as well.

[jedbrown]

Cool. Can we test isentropic vortex now or does it need some conversion? Do we have any problems in which we can show convergence under refinement?

[LeilaGhaffari]

Cool. Can we test isentropic vortex now or does it need some conversion? Do we have any problems in which we can show convergence under refinement?

Unfortunately, isentropic vortex doesn't have support for primitive variables. If we want to test the full functionality, we should either make the channel problem compute the errors (I tried but it was more complicated than switching on a flag) or, refactor the euler solver to use newtonian (I think this is the simplest way to make traveling vortex functional for our desired convergence study).

Or we can just show the convergence with isentropic vortex using conservative variables (without the conversion). This is ready!

However, I am still concerned about this issue. Could you take a look at the functions I defined/modified in fluids/misc.c (mainly FixMultiplicity_NS() and ComputeL2Error())?

[LeilaGhaffari]

@jedbrown , I think I fixed the issues you raised the other day but I think the problem I am seeing with the parallel run could be irrelevant to this PR.
When testing the blasius problem on main I get:

Serial

libCEED (main)$ build/fluids-navierstokes -options_file examples/fluids/tests-output/blasius_test.yaml -ts_max_steps 10

-- Navier-Stokes solver - libCEED + PETSc --
  MPI:
    Host Name                          : leila-ThinkPad-P53s
    Total ranks                        : 1
  Problem:
    Problem Name                       : blasius
    Stabilization                      : none
  libCEED:
    libCEED Backend                    : /cpu/self/opt/blocked
    libCEED Backend MemType            : host
  PETSc:
    Box Faces                          : 3,20,1
    DM MatType                         : aij
    DM VecType                         : standard
    Time Stepping Scheme               : implicit
  Mesh:
    Number of 1D Basis Nodes (P)       : 2
    Number of 1D Quadrature Points (Q) : 2
    Global DoFs                        : 656
    DoFs per node                      : 5
    Global 5-DoF nodes                 : 131
  Partition:                             (min,max,median,max/median)
    Global Vector 5-DoF nodes          : 131, 131, 131, 1.000000
    Local Vector 5-DoF nodes           : 168, 168, 168, 1.000000
    Ghost Interface 5-DoF nodes        : 0, 0, 0, -nan
    Ghost Interface Ranks              : 0, 0, 0, -nan
    Owned Interface 5-DoF nodes        : 0, 0, 0, -nan
    Owned Interface Ranks              : 0, 0, 0, -nan
Time taken for solution (sec): 0.224105
Time integrator CONVERGED_ITS on time step 10 with final time 2e-05``

Parallel:

libCEED (main)$ mpiexec.hydra -n 6 build/fluids-navierstokes -options_file examples/fluids/tests-output/blasius_test.yaml -ts_max_steps 10

-- Navier-Stokes solver - libCEED + PETSc --
  MPI:
    Host Name                          : leila-ThinkPad-P53s
    Total ranks                        : 6
  Problem:
    Problem Name                       : blasius
    Stabilization                      : none
  libCEED:
    libCEED Backend                    : /cpu/self/opt/blocked
    libCEED Backend MemType            : host
  PETSc:
    Box Faces                          : 3,20,1
    DM MatType                         : aij
    DM VecType                         : standard
    Time Stepping Scheme               : implicit
  Mesh:
    Number of 1D Basis Nodes (P)       : 2
    Number of 1D Quadrature Points (Q) : 2
    Global DoFs                        : 656
    DoFs per node                      : 5
    Global 5-DoF nodes                 : 131
  Partition:                             (min,max,median,max/median)
    Global Vector 5-DoF nodes          : 17, 28, 20, 1.384615
    Local Vector 5-DoF nodes           : 36, 36, 36, 1.000000
    Ghost Interface 5-DoF nodes        : 0, 8, 8, 1.000000
    Ghost Interface Ranks              : 0, 1, 1, 1.000000
    Owned Interface 5-DoF nodes        : 0, 8, 8, 1.000000
    Owned Interface Ranks              : 0, 1, 1, 1.000000
Time taken for solution (sec): 0.353972
Time integrator DIVERGED_NONLINEAR_SOLVE on time step 0 with final time 0.

Basically, the solver is diverging with parallel runs while there is no problem with the serial run. Is this normal?
(There could still be some issues with my implementation though. Could you give another look at your convenience?)

Oh, I shouldn't have squashed the commits, I though it could help. You can see your suggestions here.

[jedbrown]

Why is the nonlinear solve diverging in parallel? Is the linear solver failing? Is this solely a difference of the preconditioner getting weaker in parallel for the steady problem?

[LeilaGhaffari]

The failure is very random and it happens in the linear solver. I ran some tests on Noether with main. For example, it failed with 6/8/20 processors but not with 1-5/7/10!
Here is an example of the failed ones:

$ mpiexec -n 8 build/fluids-navierstokes -options_file examples/fluids/tests-output/blasius_test.yaml -ts_max_steps 10 -ksp_converged_reason -snes_converged_reason

-- Navier-Stokes solver - libCEED + PETSc --
  MPI:
    Host Name                          : noether
    Total ranks                        : 8
  Problem:
    Problem Name                       : blasius
    Stabilization                      : none
  libCEED:
    libCEED Backend                    : /cpu/self/opt/blocked
    libCEED Backend MemType            : host
  PETSc:
    Box Faces                          : 3,20,1
    DM MatType                         : aij
    DM VecType                         : standard
    Time Stepping Scheme               : implicit
  Mesh:
    Number of 1D Basis Nodes (P)       : 2
    Number of 1D Quadrature Points (Q) : 2
    Global DoFs                        : 656
    DoFs per node                      : 5
    Global 5-DoF nodes                 : 131
  Partition:                             (min,max,median,max/median)
    Global Vector 5-DoF nodes          : 12, 22, 16, 1.400000
    Local Vector 5-DoF nodes           : 28, 32, 28, 1.142857
    Ghost Interface 5-DoF nodes        : 0, 8, 8, 1.000000
    Ghost Interface Ranks              : 0, 1, 1, 1.000000
    Owned Interface 5-DoF nodes        : 0, 8, 8, 1.000000
    Owned Interface Ranks              : 0, 1, 1, 1.000000
    Linear solve did not converge due to DIVERGED_ITS iterations 10000
  Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
Time taken for solution (sec): 0.257922
Time integrator DIVERGED_NONLINEAR_SOLVE on time step 0 with final time 0.

(I don't see any preconditioning being set in the yaml file.)

[jedbrown]

You can use -ts_view -ksp_converged_reason -ksp_view_singularvalues to better understand the failure. Does the same failure occur on main?

[LeilaGhaffari]

Oh, yes. Those runs and the following are all on main.

libCEED (main)$ mpiexec.hydra -n 6 build/fluids-navierstokes -options_file examples/fluids/tests-output/blasius_test.y
aml -ts_max_steps 10 -ts_view -ksp_converged_reason -ksp_view_singularvalues

-- Navier-Stokes solver - libCEED + PETSc --
  MPI:
    Host Name                          : leila-ThinkPad-P53s
    Total ranks                        : 6
  Problem:
    Problem Name                       : blasius
    Stabilization                      : none
  libCEED:
    libCEED Backend                    : /cpu/self/opt/blocked
    libCEED Backend MemType            : host
  PETSc:
    Box Faces                          : 3,20,1
    DM MatType                         : aij
    DM VecType                         : standard
    Time Stepping Scheme               : implicit
  Mesh:
    Number of 1D Basis Nodes (P)       : 2
    Number of 1D Quadrature Points (Q) : 2
    Global DoFs                        : 656
    DoFs per node                      : 5
    Global 5-DoF nodes                 : 131
  Partition:                             (min,max,median,max/median)
    Global Vector 5-DoF nodes          : 17, 28, 20, 1.384615
    Local Vector 5-DoF nodes           : 36, 36, 36, 1.000000
    Ghost Interface 5-DoF nodes        : 0, 8, 8, 1.000000
    Ghost Interface Ranks              : 0, 1, 1, 1.000000
    Owned Interface 5-DoF nodes        : 0, 8, 8, 1.000000
    Owned Interface Ranks              : 0, 1, 1, 1.000000
    Linear solve did not converge due to DIVERGED_ITS iterations 10000
Iteratively computed extreme singular values: max 482.958 min 1.45981 max/min 330.837
TS Object: 6 MPI processes
  type: beuler
  maximum steps=10
  maximum time=500.
  total number of I function evaluations=1
  total number of I Jacobian evaluations=1
  total number of nonlinear solver iterations=0
  total number of linear solver iterations=10000
  total number of nonlinear solve failures=1
  total number of rejected steps=1
  using relative error tolerance of 0.0001,   using absolute error tolerance of 0.0001
  TSAdapt Object: 6 MPI processes
    type: none
  SNES Object: 6 MPI processes
    type: newtonls
    maximum iterations=50, maximum function evaluations=10000
    tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
    total number of linear solver iterations=10000
    total number of function evaluations=1
    norm schedule ALWAYS
    SNESLineSearch Object: 6 MPI processes
      type: bt
        interpolation: cubic
        alpha=1.000000e-04
      maxstep=1.000000e+08, minlambda=1.000000e-12
      tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08
      maximum iterations=40
    KSP Object: 6 MPI processes
      type: gmres
        restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
        happy breakdown tolerance 1e-30
      maximum iterations=10000, initial guess is zero
      tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
      left preconditioning
      using PRECONDITIONED norm type for convergence test
    PC Object: 6 MPI processes
      type: bjacobi
        number of blocks = 6
        Local solver information for first block is in the following KSP and PC objects on rank 0:
        Use -ksp_view ::ascii_info_detail to display information for all blocks
      KSP Object: (sub_) 1 MPI process
        type: preonly
        maximum iterations=10000, initial guess is zero
        tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
        left preconditioning
        using NONE norm type for convergence test
      PC Object: (sub_) 1 MPI process
        type: ilu
          out-of-place factorization
          0 levels of fill
          tolerance for zero pivot 2.22045e-14
          matrix ordering: natural
          factor fill ratio given 1., needed 1.
            Factored matrix follows:
              Mat Object: (sub_) 1 MPI process
                type: seqaij
                rows=88, cols=88
                package used to perform factorization: petsc
                total: nonzeros=3680, allocated nonzeros=3680
                  using I-node routines: found 26 nodes, limit used is 5
        linear system matrix = precond matrix:
        Mat Object: (sub_) 1 MPI process
          type: seqaij
          rows=88, cols=88
          total: nonzeros=3680, allocated nonzeros=0
          total number of mallocs used during MatSetValues calls=0
            using I-node routines: found 26 nodes, limit used is 5
      linear system matrix = precond matrix:
      Mat Object: Mat_0x84000004_0 6 MPI processes
        type: mpiaij
        rows=656, cols=656
        total: nonzeros=37920, allocated nonzeros=0
        total number of mallocs used during MatSetValues calls=0
          using I-node (on process 0) routines: found 26 nodes, limit used is 5
Time taken for solution (sec): 0.333688
Time integrator DIVERGED_NONLINEAR_SOLVE on time step 0 with final time 0.

[jeremylt]

Closing in anticipation of moving Honee to a new repo

this seems like a good branch to recreate on the other side of the move! We should make an issue when we stand up the new repo @jrwrigh