Picard solver as alternative to Newton for iNS #397

richardschu · 2023-03-28T14:53:00Z

Implemented a Picard solver as an alternative to Newton for the incompressible flow problems, since it might be interesting to do only a fixed single nonlinear iteration in time-dependent (not quasi-static) problems. Within a Picard solver, this is a sensible linearization in time, whereas for Newton, it still is a single nonlinear step (possibly equally valid). The Poiseuille flow application was used for testing and changes to it will be removed before merging. I tested L2 errors in the velocity and pressure, linear-, and nonlinear iteration counts for: PressureInflow/ParabolicInflow, always with symmetry boundary, Picard or Newton, 'BDFCoupledSolution' or 'BDFPressureCorrection', 'DivergenceFormulation' or 'ConvectiveFormulation' for the convective term.
-All 32 cases give comparable errors (~1e-10 or less), Newton needs less nonlinear iterations, the linear solver iteration counts sometimes slightly increase for Picard, but the convergence criteria are identical in both cases (and the systems are vastly different, so that might in fact not be too easy to compare)

The structure module uses the nonlinear solver as well, but the operators are not implemented (assertions placed). I think Picard does not make too much sense for structure problems, at least it is unusual.

The filename newton_solver.h/_data.h is kept for comparison and will be changed.

A lot of the changes are only due to default 'NonlinearSolver::SolverType::Newton' (which I did not want to hide) and renaming in all the modules.

I changed the tests in the Poiseuille flow application - in a follow-up PR, I will try to test more combinations. We should think about how we want to unit-test!

richardschu · 2023-03-28T14:55:52Z

applications/fluid_structure_interaction/bending_wall/application.h

@@ -249,7 +249,8 @@ class Application : public FluidFSI::ApplicationBase<dim, Number>
    // momentum step

    // Newton solver
-    param.newton_solver_data_momentum = Newton::SolverData(100, ABS_TOL, REL_TOL);
+    param.nonlinear_solver_data_momentum =
+      NonlinearSolver::SolverData(100, ABS_TOL, REL_TOL, NonlinearSolver::SolverType::Newton);


All applications define the nonlinear solver with a type to make the available choice obvious.
Alternatively, one can define NonlinearSolver::SolverType::Newton as default and omit it in all the applications.

richardschu · 2023-03-28T14:57:25Z

applications/incompressible_navier_stokes/poiseuille/application.h

@@ -196,6 +196,18 @@ class Application : public ApplicationBase<dim, Number>
                        apply_symmetry_bc,
                        "Apply symmetry boundary condition.",
                        dealii::Patterns::Bool());
+      prm.add_parameter("UseNewtonElsePicard",


these will be removed, just so you see how i tested.

richardschu · 2023-03-28T14:58:53Z

applications/incompressible_navier_stokes/poiseuille/application.h

@@ -264,12 +284,12 @@ class Application : public ApplicationBase<dim, Number>

    // pressure Poisson equation
    this->param.solver_pressure_poisson         = SolverPressurePoisson::CG;
-    this->param.solver_data_pressure_poisson    = SolverData(1000, 1.e-20, 1.e-6, 100);
+    this->param.solver_data_pressure_poisson    = SolverData(1000, 1.e-12, 1.e-4, 100);


I had to increase the tolerance, since Picard has no quadratic convergence without acceleration.
These values are not too bad, though.

why didn't you increase the maximum number of iterations instead?

Because testing is slower using more iterations and the problem we are looking at is quasi-stationary.
I the solution does not change, 1e-6 relative is tough and 1e-20 absolute is not a lot either.
If I lower the linear solver tolerances again, I can increase the tolerances in the nonlinear solver I suppose, increasing the max. number of iterations.

richardschu · 2023-03-28T14:59:56Z

applications/incompressible_navier_stokes/poiseuille/application.h


  FormulationViscousTerm const formulation_viscous_term =
    FormulationViscousTerm::LaplaceFormulation;

-  double const max_velocity = 1.0;
+  double const max_velocity = 1.0e-1;


this was lowered for the 'PressureInflow' tests. Otherwise we actually need backflow stabilization, which messes up the error we compute.

richardschu · 2023-03-28T15:00:27Z

applications/incompressible_navier_stokes/poiseuille/input.json

@@ -7,16 +7,19 @@
    "SpatialResolution": {
        "DegreeMin": "2",
        "DegreeMax": "2",
-        "RefineSpaceMin": "0",
-        "RefineSpaceMax": "0"
+        "RefineSpaceMin": "1",


Changes to the input.json here will be removed completely - just to show how I tested.

richardschu · 2023-03-28T15:18:13Z

include/exadg/solvers_and_preconditioners/nonlinear_solvers/newton_solver.h

@@ -19,18 +19,18 @@
 *  ______________________________________________________________________


TODO: rename the files.

richardschu · 2023-03-28T15:18:47Z

include/exadg/solvers_and_preconditioners/nonlinear_solvers/newton_solver.h

-    increment.reinit(solution);
-    temporary.reinit(solution);
+    VectorType residual_rhs, nonlinear_solver_output, temporary_solution;
+    residual_rhs.reinit(solution);


holds residual/rhs. I did not want to introduce yet another vector.

richardschu · 2023-03-28T15:21:24Z

include/exadg/solvers_and_preconditioners/nonlinear_solvers/newton_solver.h

-                  "Maximum number of iterations exceeded!"));
+                  "Nonlinear solver failed to solve nonlinear problem to given tolerance. "
+                  "Maximum number of iterations (" +
+                  std::to_string(solver_data.max_iter) + ") exceeded!"));


If we want to use linearized solvers, we have to give the nonlinear solvers additional parameters: 'enable_damping' (and more parameters potentially, these are hard-coded above), and some switch to activate the Assertion here, such that doing 5 nonlinear iters is ok, no matter the reduction in the residual.

richardschu · 2023-03-28T15:23:17Z

include/exadg/structure/spatial_discretization/operator.cpp

@@ -304,6 +304,21 @@ Operator<dim, Number>::setup_operators()
  operator_data.material_descriptor = material_descriptor;
  operator_data.unsteady            = (param.problem_type == ProblemType::Unsteady);
  operator_data.density             = param.density;
+


The structure module has no Picard yet. The nonlinear solver requires the according data structure and the rhs for Picard nonetheless. Both are asserted for 'NonlinearSolver::SolverType::Picard'.

richardschu · 2023-03-28T15:23:44Z

include/exadg/structure/spatial_discretization/operator.h

+  {
+    AssertThrow(
+      false, dealii::ExcMessage("Nonlinear Picard solver not implemented for structure module."));
+  }


Picard for elasticity might be useless.

…_correction.cpp

nfehn · 2023-03-28T17:01:22Z

include/exadg/solvers_and_preconditioners/nonlinear_solvers/newton_solver.h

+      // reset iterate and rhs of nonlinear solver
+      if(solver_data.nonlinear_solver_type == SolverType::Newton)
+      {
+        // zero increment as initial guess for Newton
+        nonlinear_solver_output = 0.0;
+
+        // overwrite residual with right-hand side for Newton solver
+        // multiply by -1.0 since the linearized problem is
+        // "linear_operator * increment = - residual"
+        residual_rhs *= -1.0;
+      }
+      else if(solver_data.nonlinear_solver_type == SolverType::Picard)
+      {
+        // last solution as initial guess for Picard solver
+        nonlinear_solver_output = solution;


I am not convinced that this "if-else-logic" is how we should implement a Picard solver. Generally, one should be able to introduce a Picard solver without touching every application. The current realization does not seem to be modular in this regard.

One does not have to touch all applications, Newton is default and is used if not specifically requested otherwise - I just did it to specifically show that it can be requested. Otherwise it is somehow 'hidden' for some hypothetical user.
Of course I can go the route of introducing a base class and derived ones ... but I am not sure if that helps. If we had some completely different scheme being implemented by someone (I am not even familiar with more alternative schemes), then it would make sense to come up with such a design. Keep in mind, though, that such alternative solution schemes might need entirely different operators/steps and whatnot. For Picard and Newton this currently fits - if we want to extend later, we can find a better design that fits the purpose.

If am of a different opinion, because making changes will become more and more difficult. The present solution is in my opinion a quick solution but not a clean design.

So we should/could cook up a design for nonlinear solvers in general in base/derived classes. From my perspective, I can only think of solvers similar to Newton & Picard, since I am not familiar with other methods. I assume that other methods might even need other functionality than what we have in, say, 'convective_operator.h/.cpp'. I can generalize the nonlinear solver, but still targeting Newton or Picard.

nfehn · 2023-03-28T17:10:26Z

I understand the description of this PR as introducing a new numerical method whose benefit has not been demonstrated so far? Is the idea to have this type of solver in a kind of "experimental mode" in ExaDG? Against this background, touching 70 files is immense and I wonder whether this is the way to go. For the future, it would be better to first discuss and agree on a basic design/realization of a Picard solver (touching only a single file) before making changes with 1500 new lines of code in 70 files.

nfehn · 2023-03-28T17:15:07Z

include/exadg/structure/user_interface/enum_types.h

+enum class LinearizationType
+{
+  Undefined,
+  Newton
+};


What is the meaning of this enum if it only allows one useful value?

We might either i) use the Undefined as default, or ii) someone wants to define a new nonlinear solver. Not the way to go?

richardschu · 2023-03-28T17:16:12Z

It is meant as an alternative, if I come up with some example, that is even more changes. I wanted to have a base version first. The tests I did even show that it is inferior to Newton, but maybe in some cases it is not? That should be ok, shouldn't it?
I thought that when I did not do all the changes (renaming, default in applications = at least 50% of the files changed), we would be discussing them.

nfehn · 2023-03-28T20:08:56Z

Let me explain from a different perspective why the current design is problematic:

Because of the call nonlinear_operator.evaluate_rhs_picard(residual_rhs);, every nonlinear operator has to implement now a function evaluate_rhs_picard(), even though some module might never want to use Picard!
Since you plug the nonlinear solver type into the SolverData, every high-level module using any type of nonlinear solver now depends on this parameter, even though a certain module might never want to use Picard!
the circumstance that you renamed variables in the Newton solver indicates that combining the two cases in a single implementation in the current form is not intuitive and hinders understanding of the code.

In my opinion, this is not a modular design, because introducing a new feature for one module implies making changes to all other modules, even though they might never want to use this feature. It is not in line with the design that we want to follow in ExaDG in my opinion. In this respect, I would argue that we want to follow the library idea of deal.II also in ExaDG. Take as an example the implementation (or wrappers) for linear solvers. According to your design, implementing a new linear solver for non-symmetric problems would imply making changes to the Poisson problem with symmetric positive definite linear operator!? (Note that "making changes" could also mean "recompiling the code", since the nonlinear solver type parameter is included in the nonlinear solver implementation, which is then included by all modules).

If I remember correctly, this topic is also discussed under the name "dependency inversion principle" in software design ("high-level modules should depend on abstractions", abstraction meaning something like "interface", e.g. realized via a base class with inheritance and dynamic polymorphism). I do not want to imply that inheritance is the only option. The important aspect - according to my understanding - is to realize the very general goal in software development of minimizing dependencies. So far, I only thought of Newton as a nonlinear solver. So the argument that we need to do the changes only once and for all and never touch it in the future is not convincing in my opinion. By such an argument actually, other project copied the Newton solver implementation from ExaDG, arguing that a Newton solver is written once and then never changed. (For clarity, I would not change Newton because of a new Picard implementation, but my point is that there are multiple possibilities to extend the implementation of an existing Newton solver itself. So my point is one should probably avoid to argue with "never" in software development ... In most cases when I had taken this perspective in the past, I have proved myself wrong in this regard.)

I hope these explanations help, but we can of course also discuss this topic in person.

richardschu · 2023-03-29T07:22:50Z

Hi Niklas,
thank you for your time and effort checking out the PR. I see the points you are making, sounds reasonable and I even experienced that when I had to implement a evaluate_rhs_picard in the structure solver. The current state of the PR is pure functionality, so that we have a common starting ground and I could get familiar with that part of the code.
I will try to come up with a better design and keep you posted!

richardschu added 4 commits March 27, 2023 19:28

picard solver added, tested with poiseuille flow

cc41d08

renaming and typo fix

1e71944

renaming and merging in convective_operator

6d6e6d1

indentation

88edd7f

richardschu commented Mar 28, 2023

View reviewed changes

richardschu added 3 commits March 28, 2023 17:27

removed the commented dst.add(1.0, *rhs_vector); in operator_pressure…

5dcda26

…_correction.cpp

adapted tests for Poiseuille

1ccb47c

indentation

feba993

nfehn reviewed Mar 28, 2023

View reviewed changes

nfehn marked this pull request as draft June 2, 2023 13:08

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Picard solver as alternative to Newton for iNS #397

Picard solver as alternative to Newton for iNS #397

richardschu commented Mar 28, 2023 •

edited

richardschu Mar 28, 2023

richardschu Mar 28, 2023

richardschu Mar 28, 2023

nfehn Mar 28, 2023

richardschu Mar 28, 2023

richardschu Mar 28, 2023

richardschu Mar 28, 2023

richardschu Mar 28, 2023

richardschu Mar 28, 2023

richardschu Mar 28, 2023

richardschu Mar 28, 2023

richardschu Mar 28, 2023

nfehn Mar 28, 2023

richardschu Mar 28, 2023

nfehn Mar 28, 2023

richardschu Mar 28, 2023

nfehn commented Mar 28, 2023

nfehn Mar 28, 2023

richardschu Mar 28, 2023

richardschu commented Mar 28, 2023

nfehn commented Mar 28, 2023

richardschu commented Mar 29, 2023

		@@ -19,18 +19,18 @@
		* ______________________________________________________________________

Picard solver as alternative to Newton for iNS #397

Are you sure you want to change the base?

Picard solver as alternative to Newton for iNS #397

Conversation

richardschu commented Mar 28, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nfehn commented Mar 28, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

richardschu commented Mar 28, 2023

nfehn commented Mar 28, 2023

richardschu commented Mar 29, 2023

richardschu commented Mar 28, 2023 •

edited