Possible causes of "Loss of Precision" error from Aztec00 in turbulent 2D flow over a cylinder?

[ShadowEngineer]

I'm attempting to run a simple 2D incompressible flow over a cylinder (similar to the example) simulation for a prolonged period of time (at least 3 minutes of sim time) in very turbulent conditions at 50000 Reynolds number with initial velocity 1 and kinematic viscosity 2e-05 with:

1. with gls_navier_stokes_2d solver
2. my own mesh built with gmsh
3. adaptive time stepping with a target CFL of 2
4. mesh adaption with max 1_000_000 mesh elements
5. amg linear solver with minimum residual 1e-12
6. inexact newton non-linear solver

After many re-runs and attempts with various results, I keep getting the same error:

AztecOO::Iterate error code -3: loss of precision

After diving through the Trilinos source code, I found this definition in az_solve.c

  status:          On output, indicates termination status:
                    0:  terminated normally.
                   -1:  maximum number of iterations taken without achieving
                        convergence.
                   -2:  Breakdown. The algorithm can not proceed due to
                        numerical difficulties (usually a divide by zero).
                   -3:  Internal residual differs from the computed residual due
                        to a significant loss of precision.

I've attempted to change many variables in my parameters file, including:

1. lowering minimum residual in linear solver (it used to be 1e-8, and is now 1e-12)
2. increasing max mesh elements in refinement (used to be 100k, now is 1 mil)
3. enabling mesh refinement in the first place (i started with it off)
4. enabling adaptive time stepping in the first place (i also started with it off)

All of these only made the simulation survive for longer, not fix the issue completely. I've tried looking through lethe documentation and deal.II and trilinos documentation to better understand what might be causing this error, but have not found anything of great use.

So my questions are:

1. Is this a known issue in lethe?
2. Is it possible to fix this problem completely?
3. If so, what can I do to fix it?
4. If no, how can I mitigate it more, or is what I've done the best I can do?
5. Is this the best place to ask about this?
6. Could this be a bug?

A video of the simulation has been made, for demonstration purposes. At the 0:10 second mark, notice the large growth of >2 m/s velocity after a vortex tries to leave the Outlet right boundary, causing a rapid growth of turbulent kinetic energy in the system and a subsequent "loss of precision" error causing the simulation to end.

My .geo file for generating the mesh and the .prm for the simulation are attached:

// CONFIGURATION
PIPE_WIDTH = 50;
PIPE_HEIGHT = 30;
CYLINDER_X = 10;
CYLINDER_Y = 15;
CYLINDER_RADIUS = 0.5;
LINE_PROGRESSION_RANGE_1 = 15; // 30
LINE_PROGRESSION_RANGE_2 = 15; // 30
LINE_PROGRESSION_RANGE_3 = 80; // 150
TARGET_MESH_SIZE = 0.20;

// variables
y_coord_upper = 0;
y_coord_lower = PIPE_HEIGHT;
x_coord_lower = 0;
x_coord_upper = PIPE_WIDTH;

root_cylinder_radius = 1.414213562 / 2. * CYLINDER_RADIUS;
cylinder_radius_long = 10 * CYLINDER_RADIUS;
root_cylinder_radius_long = 1.414213562 / 2. * cylinder_radius_long;

cylinder_border_left = CYLINDER_X - root_cylinder_radius_long;
cylinder_border_right = CYLINDER_X + root_cylinder_radius_long;
cylinder_border_top = CYLINDER_Y - root_cylinder_radius_long;
cylinder_border_bottom = CYLINDER_Y + root_cylinder_radius_long;

// top and bottom wall points
Point(0) = {x_coord_lower,  y_coord_upper,  0,  TARGET_MESH_SIZE};
Point(1) = {x_coord_lower,  y_coord_lower,  0,  TARGET_MESH_SIZE};
Point(2) = {cylinder_border_left,  y_coord_upper,  0,  TARGET_MESH_SIZE};
Point(3) = {cylinder_border_left,  y_coord_lower,  0,  TARGET_MESH_SIZE};
Point(4) = {cylinder_border_right,  y_coord_upper,  0,  TARGET_MESH_SIZE};
Point(5) = {cylinder_border_right,  y_coord_lower,  0,  TARGET_MESH_SIZE};
Point(6) = {x_coord_upper,  y_coord_upper,  0,  TARGET_MESH_SIZE};
Point(7) = {x_coord_upper,  y_coord_lower,  0,  TARGET_MESH_SIZE};

// left and right wall intermediate points
Point(11) = {x_coord_lower,  CYLINDER_Y - root_cylinder_radius_long,  0,  TARGET_MESH_SIZE};
Point(12) = {x_coord_lower,  CYLINDER_Y + root_cylinder_radius_long,  0,  TARGET_MESH_SIZE};
Point(13) = {x_coord_upper,  CYLINDER_Y - root_cylinder_radius_long,  0,  TARGET_MESH_SIZE};
Point(14) = {x_coord_upper,  CYLINDER_Y + root_cylinder_radius_long,  0,  TARGET_MESH_SIZE};

// cylinder points
Point(20) = {CYLINDER_X,  CYLINDER_Y,  0,  TARGET_MESH_SIZE};
Point(21) = {CYLINDER_X - root_cylinder_radius,  CYLINDER_Y - root_cylinder_radius,  0,  TARGET_MESH_SIZE};
Point(22) = {CYLINDER_X + root_cylinder_radius,  CYLINDER_Y - root_cylinder_radius,  0,  TARGET_MESH_SIZE};
Point(23) = {CYLINDER_X + root_cylinder_radius,  CYLINDER_Y + root_cylinder_radius,  0,  TARGET_MESH_SIZE};
Point(24) = {CYLINDER_X - root_cylinder_radius,  CYLINDER_Y + root_cylinder_radius,  0,  TARGET_MESH_SIZE};

// cylinder border corners
Point(31) = {cylinder_border_left,  cylinder_border_top,  0,  TARGET_MESH_SIZE};
Point(32) = {cylinder_border_right,  cylinder_border_top,  0,  TARGET_MESH_SIZE};
Point(33) = {cylinder_border_right,  cylinder_border_bottom,  0,  TARGET_MESH_SIZE};
Point(34) = {cylinder_border_left,  cylinder_border_bottom,  0,  TARGET_MESH_SIZE};

// lines between all the points
Line(60) = {0, 11};
Line(61) = {11, 12};
Line(63) = {12, 1};
Line(64) = {6, 13}; // right wall top segment
Line(65) = {13, 14}; // right wall middle segment
Line(66) = {14, 7}; // right wall bottom segment

Line(1) = {0, 2}; // left segment of top wall
Line(2) = {1, 3}; // left segmnet of bottom wall
Line(3) = {2, 4}; // middle segment of top wall
Line(4) = {3, 5}; // middle segment of bottom wall
Line(5) = {4, 6}; // right segment of top wall
Line(6) = {5, 7}; // right segment of bottom wall

Line(10) = {11, 12};
Line(11) = {11, 31};
Line(12) = {12, 34};
Line(13) = {31, 34};

Line(20) = {32, 33};
Line(21) = {32, 13};
Line(22) = {13, 14};
Line(23) = {33, 14};

Line(30) = {2, 31};
Line(31) = {4, 32};
Line(32) = {31, 32};
Line(33) = {34, 3};
Line(34) = {34, 33};
Line(35) = {33, 5};

Circle(40) = {24, 20, 23};
Circle(41) = {23, 20, 22};
Circle(42) = {22, 20, 21};
Circle(43) = {21, 20, 24};

Line(50) = {24, 34};
Line(51) = {23, 33};
Line(52) = {22, 32};
Line(53) = {21, 31};

//Theta of the circle
Transfinite Line {3, 22, 4, 10, 34, 20, 32, 13, 40, 41, 42, 43} = Ceil(LINE_PROGRESSION_RANGE_2) Using Progression 1.0;

// Radial direction
Transfinite Line {50, 51, 52, 53} = Ceil(LINE_PROGRESSION_RANGE_1) Using Progression 1.05;
Transfinite Line {5, 21, 23, 6} = Ceil(LINE_PROGRESSION_RANGE_3) Using Progression 1.00;

Line Loop(41) = {-50, 40, 51, -34};
Line Loop(42) = {-13, -53, 43, 50};
Line Loop(43) = {53, 32, -52, 42};
Line Loop(44) = {52, 20, -51, 41};

Line Loop (50) = {-60, 1, 30, -11};
Line Loop (51) = {-10, 11, 13, -12};
Line Loop (52) = {-63, 12, 33, -2};
Line Loop (53) = {-30, 3, 31, -32};
Line Loop (54) = {-33, 34, 35, -4};
Line Loop (55) = {-31, 5, 64, -21};
Line Loop (56) = {-20, 21, 22, -23};
Line Loop (57) = {-35, 23, 66, -6};

Plane Surface(41) = {41};
Plane Surface(42) = {42};
Plane Surface(43) = {43};
Plane Surface(44) = {44};
Plane Surface(50) = {50};
Plane Surface(51) = {51};
Plane Surface(52) = {52};
Plane Surface(53) = {53};
Plane Surface(54) = {54};
Plane Surface(55) = {55};
Plane Surface(56) = {56};
Plane Surface(57) = {57};
Transfinite Surface {41:44, 50:57};
Recombine Surface {41:44, 50:57};

Physical Surface(1) = {41:44, 50:57};
Physical Line(0) = {40, 41, 42, 43}; // cylinder
Physical Line(1) = {60, 10, 63}; // left side,  inlet
Physical Line(2) = {1, 3, 5, 2, 4, 6}; // top and bottom walls
Physical Line(3) = {64, 22, 66}; // right side, outlet

Mesh 1;
Mesh 2;
Save "pipe-with-cylinder.msh";

subsection mesh
    set type = gmsh
    set file name = pipe-with-cylinder.msh
end

subsection simulation control
    set method = bdf2

    set output name = output
    set output frequency = 40
    set output path = ./Output/

    set time end = 120
    set time step = 0.01
    set subdivision = 1

    # adaptive time stepping
    set adapt = true
    set max cfl = 2
    # set max time step = 1 # as a safeguard so the simulation doesn't go crazy
    set adaptative time step scaling = 1.1
end

subsection post-processing
    set smoothed output fields = true
end

# dynamically subdivide or recombine mesh elements based on how much they are contributing to error.
# basically, the areas where there's too much detail to fit into one cell (too high error),
# mesh adaptation will divide that singular element into multiple

subsection mesh adaptation
    set type = kelly
    set variable = velocity
    set frequency = 1
    set min refinement level = 0
    set max refinement level = 3
    set max number elements = 10000000
    set fraction coarsening = 0.05
    set fraction refinement = 0.1
    set fraction type = number
    set initial refinement steps = 0
end

# simply making every node (element, I think???) originally start with a velocity of 1 m/s moving right
subsection initial conditions
    set type = nodal
    subsection uvwp
        set Function expression = 1; 0; 0
    end
end

subsection boundary conditions
    set number = 4
    set time dependent = false

    # circle (the cylinder in the way of the fluid flow)
    subsection bc 0
        set type = noslip
    end

    # left side (fluid entering the pipe)
    subsection bc 1
        set type = function

        # constantly "moving water into the system" at a steady velocity of 1 m/s right
        subsection u
            set Function expression = 1
        end

        subsection v
            set Function expression = 0
        end

        subsection w
            set Function expression = 0
        end
    end

    # top and bottom walls
    subsection bc 2
        set type = noslip
    end

    # right wall
    subsection bc 3
        set type = outlet
    end
end

subsection physical properties
    set number of fluids = 1
    subsection fluid 0 
        set kinematic viscosity = 2e-05 # kinematic viscosity of 2e-05 m^2/s
    end
end

subsection FEM
    # linear interpolation for velocity
    set velocity order = 1

    # linear interpolation for pressure too
    set pressure order = 1
end

subsection manifolds
    set number = 1
    subsection manifold 0
        set id   = 0
        set type = spherical
        set arg1 = 8
        set arg2 = 8
    end
end

# got these from the pipe-with-cylinder.prm example file from lethe
subsection non-linear solver
    set verbosity      = verbose
    set solver         = inexact_newton
    set tolerance      = 1e-6
    set max iterations = 10
end
subsection linear solver
    set verbosity                                 = verbose
    set method                                    = amg
    set relative residual                         = 1e-4
    set minimum residual                          = 1e-12
    set amg preconditioner ilu fill               = 0
    set amg preconditioner ilu absolute tolerance = 1e-12
    set amg preconditioner ilu relative tolerance = 1.00
end

subsection forces
    set calculate force = true
    set calculate torque = false
    set calculation frequency = 1
    set output frequency = 100
    set output precision = 10
    set force name = force
    set verbosity = quiet
end

subsection timer
    set type = iteration
end

[blaisb]

Dear @ShadowEngineer Sorry for the delay in getting back at you. It was a good idea to open an issue.

The solver you are reporting (gls_navier_stokes_2d) is very dated, we have changed the naming scheme of the solvers around 8 months ago. That solver is now lethe-fluid. Still the underlying numerics have not changed significantly.

The issue you are getting (as seen in your movie) that is your simulation crashes as soon as you have a vortex hitting the outlet. This is a well-established issue with "do-nothing" boundary conditions in FEM because they impose a zero traction at the outlet boundary condition. In essence, this is a blend between imposing a zero pressure value and a zero gradient of the normal velocity. The issue is that when you have backflow (aka a vortex that touches the wall) you have re-entry of fluid and the results become essentially garbage since it is unclear with what information the fluid comes back in.

However, we actually have a good way to mitigate this. The best approach is to use a the outlet boundary condition provided in Lethe: https://chaos-polymtl.github.io/lethe/documentation/parameters/cfd/boundary_conditions_cfd.html

In the outlet boundary, there is a coefficient "beta" that you can set that prevent backflow. The idea is to set it at an increasingly high value (aka 10, 1000 or 1e4) to prevent this type of backflow while ensuring the simulation remains stable.

To achieve this, you would add the following information

    subsection bc 3
        set type = outlet
        set beta = 100
    end

They are also other alternative that can work quite well:

• Extending the length of the domain behind the cylinder
• Restricting the domain behind the cylinder after a long length to ensure that there is some sort of "venturi-effect" and that vortical structures are killed.

Yet the first thing I would try is to play with the outlet boundary and increase the beta

Feel free to reach out if you have more questions in this issue. I'll be faster to answer next time :)

[LuckaBarbeauNRC]

@blaisb I also recently got as similar error. However, in my case, it was due to the smoothing post-processor having trouble converging with very messy q_criterion results. This could also be the case here since the smoothing post-processor is activated. Maybe adding the possibility to tune the linear solver for the post-processor or adding some fail-safe mecanisme when this happen due to the smoothing could also help.