search

trilinos / Trilinos

Primary repository for the Trilinos Project
https://trilinos.org/
Other
1.19k stars 565 forks source link

Using NOX Matrix-Free Solver with User-provided Preconditioning Matrix #1991

Closed guoluanjing closed 6 years ago

guoluanjing commented 6 years ago

@trilinos/<nox/epetra>

Expectations

Applying right preconditioning (with user-provided matrix P) to the current linearized equation system Ju = -F as follows: J P^-1 Pu = -F, and use Matrix Free Newton Krylov solver to solve it in three steps: 1) linear solver to solve P^-1 v = y, or equivalently Py = v, for y; 2) approximate matrix vector product J P^-1 v = [ F(u + epsilon P^-1 v) - F(u)] / epsilon; and 3) solve the preconditioned system as a linear combination of krylov vectors: u = u0 + summation(alpha (JP^-1)^i v).

Current Behavior

We have used the nox::epetra::MatrixFree solver to solve our nonlinear equation system. But it's been very slow. So we wanted to use parts of the Jacobian matrix as the preconditioning matrix. We changed the nox::epetra::LinearSystemAztecOO constructor to NOX::Epetra::LinearSystemAztecOO::LinearSystemAztecOO(Teuchos::ParameterList & printingParams, Teuchos::ParameterList & linearSolverParams, const Teuchos::RCP< NOX::Epetra::Interface::Required > & i, const Teuchos::RCP< NOX::Epetra::Interface::Preconditioner > & iPrec, const Teuchos::RCP< Epetra_Operator > & M, const NOX::Epetra::Vector & cloneVector, const Teuchos::RCP< NOX::Epetra::Scaling > scalingObject = Teuchos::null) .

We provided the interface class that defines computePreconditioner(const Epetra_Vector&, Epetra_Operator & M, Teuchos::ParameterList* = 0)as iPrec, and the Epetra_CrsMatrix object M as the M. But the run went at the same speed, and in the output within each Nonlinear Solver Step we see the following message:

    ***** Problem: NOX::Matrix-Free
    ***** Preconditioned GMRES solution
    ***** No preconditioning
    ***** No scaling
    *******************************************************.

It seems like no preconditioning was applied whatsoever.

Motivation and Context

We are doing multi-phase fluid flow in fractured medium with applications to oil and gas reservoir simulations. Due to the phase behavior and dependencies of flow properties as well as rock properties on state variables, building an explicit Jacobian matrix is not only a daunting task but also error-prone. Matrix Free Newton Krylov nonlinear solver is very appealing. We are very happy that the implementation of the nox::epetra::MatrixFree solver without providing any Jacobian was enough to run the test problem. But the efficiency is what constrains us. It took a whole week to run a simulation for 1200 days. So we wrote part of the Jacobian to use as the preconditioning matrix. I've searched online for an example of how this is done with the MatrixFree solver in nox, but I didn't find any. Multi-phase problems exist in many application areas. When we are done with this, we could provide an example for other people to refer to.

Definition of Done

We would like to be pointed in the right direction. We have the matrix constructed, how could we use it now with the Matrix-Free solver? What are the knobs we need to turn to turn on this functionality? Or if you have an example of using user-provided preconditioning matrix with MatrixFree solver, please share it with us.

Possible Solution

In order to demonstrate how it is done, you could use the analytical Jacobian you build in the nox/test/epetra/1Dfem example as a preconditioning matrix and use it with the MatrixFree solver. This way nonlinear solve should converge with 1 step, and linear solve should converge with the same number of steps as if the conventional Newton solver was used.

Steps to Reproduce

Your Environment

Related Issues

Additional Information

aprokop commented 6 years ago

@trilinos/nox @rppawlo

rppawlo commented 6 years ago

When you supply the preconditioner operator, do you plan to do the inverse operation yourself or do you want to supply just the approximate Jacobian matrix and then use an internal preconditioner such as ifpack of ml for the the inverse? Either use case should work fine. The latter does not have an example for the epetra adapters, but does for thyra adapters. I will add one today to the 1dfem example in epetra. Could you also dump the parameter list values at the end of the run and send the output? Would like to see what defaults are being set and used.

guoluanjing commented 6 years ago

@rppawlo We would like to only supply the approximate Jacobian matrix and use an internal preconditioner. An example would be great.

Here is the parameter list:

-- Parameters Passed to Nonlinear Solver --

 Nonlinear Solver = Line Search Based
 Printing -> 
  MyPID = 0
  Output Precision = 3
  Output Processor = 0
  Output Information = 4351
  Output Stream = Teuchos::RCP<std::__1::basic_ostream<char, std::__1::char_traits<char> > >{ptr=0x7ff6ea6183c0,node=0x7ff6ea57ac60,strong_count=11,weak_count=0}
  Error Stream = Teuchos::RCP<std::__1::basic_ostream<char, std::__1::char_traits<char> > >{ptr=0x7ff6ea6183c0,node=0x7ff6ea57acd0,strong_count=6,weak_count=0}
 Direction -> 
  Method = Newton
  Newton -> 
   Forcing Term Method = Constant
   Rescue Bad Newton Solve = 1
   Linear Solver -> 
    Max Iterations = 500
    Tolerance = 1e-08
    Write Linear System = 0
    Use Preconditioner as Solver = 0   [default]
    Write Linear System File Prefix = NOX_LinSys   [default]
    Output -> 
     Total Number of Linear Iterations = 2258823   [unused]
     Number of Linear Iterations = 136   [unused]
     Achieved Tolerance = 0.178456
 Linear Solver -> 
  Jacobian Operator = Matrix-Free
  Zero Initial Guess = 1
  Aztec Solver = GMRES
  Output Frequency = 1
  Compute Scaling Manually = 0
  Preconditioner = None   [default]
  Output Solver Details = 1   [default]
  Throw Error on Prec Failure = 1   [default]
  RCM Reordering = Disabled   [default]
  Orthogonalization = Classical   [default]
  Size of Krylov Subspace = 300   [default]
  Convergence Test = r0   [default]
  Preconditioner Reuse Policy = Rebuild   [default]
  Max Age Of Prec = 1   [default]
 Line Search -> 
  Method = Backtrack
  Backtrack -> 
   Minimum Step = 1e-12   [default]
   Default Step = 1   [default]
   Recovery Step = 1   [default]
   Max Iters = 100   [default]
   Reduction Factor = 0.5   [default]
 Solver Options -> 
  Status Test Check Type = Minimal   [default]
 Output -> 
  Nonlinear Iterations = 4   [unused]
  2-Norm of Residual = 4.25848e-12   [unused]

-- Status Test Results -- ...........OR Combination -> ...........Number of Iterations = 0 < 20 ...........Finite Number Check (Two-Norm F) = Finite ...........AND Combination -> **...........WRMS-Norm = 1.000e+12 < 1 (Min Step Size: 1.000e+00 >= 1) (Max Lin Solv Tol: 1.785e-01 < 0.5) ??...........F-Norm = 0.000e+00 < 1.000e-08 (Length-Scaled Two-Norm, Absolute Tolerance)

rppawlo commented 6 years ago

From the parameter list it looks like the preconditioner is disabled. Your "Preconditioner" is set to "None" in the linear solver sublist. Try setting that to "New Ifpack". There is an example of this use case in the file:

Trilinos/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C

except that the preconditioner is set to "None" in this example too. To enable the preconditioner in this example, make this change below:

[rppawlo@gge 1Dfem]$ git diff 1DfemLeanMatrixFree.C
diff --git a/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C b/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C
index 6b20ddc..12353a8 100644
--- a/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C
+++ b/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C
@@ -197,9 +197,9 @@ int main(int argc, char *argv[])
     lsParams.set("Tolerance", 1e-4);

     // Various Preconditioner options
-    lsParams.set("Preconditioner", "None");
+    //lsParams.set("Preconditioner", "None");
     //lsParams.set("Preconditioner", "AztecOO");
-    //lsParams.set("Preconditioner", "New Ifpack");
+    lsParams.set("Preconditioner", "New Ifpack");
     lsParams.set("Preconditioner Reuse Policy", "Reuse");
     //lsParams.set("Preconditioner Reuse Policy", "Recompute");
     //lsParams.set("Preconditioner Reuse Policy", "Rebuild");
rppawlo commented 6 years ago

Resending the diff - didn't protect it properly in the issue.

diff --git a/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C b/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C
index 6b20ddc..12353a8 100644
--- a/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C
+++ b/packages/nox/test/epetra/1Dfem/1DfemLeanMatrixFree.C
@@ -197,9 +197,9 @@ int main(int argc, char *argv[])
     lsParams.set("Tolerance", 1e-4);

     // Various Preconditioner options
-    lsParams.set("Preconditioner", "None");
+    //lsParams.set("Preconditioner", "None");
     //lsParams.set("Preconditioner", "AztecOO");
-    //lsParams.set("Preconditioner", "New Ifpack");
+    lsParams.set("Preconditioner", "New Ifpack");
     lsParams.set("Preconditioner Reuse Policy", "Reuse");
     //lsParams.set("Preconditioner Reuse Policy", "Recompute");
     //lsParams.set("Preconditioner Reuse Policy", "Rebuild");
guoluanjing commented 6 years ago

Hi Roger @rppawlo, I tried to set the preconditioner to be "New Ifpack", the simulation that was converging fine before has now failed to converge. I'm sending you the current output.

Should I reset the parameters in the preconditioner? And if yes, how would I go about doing that?

Here are the parameter list values:

-- Parameters Passed to Nonlinear Solver --

 Nonlinear Solver = Line Search Based
 Printing -> 
  MyPID = 0
  Output Precision = 3
  Output Processor = 0
  Output Information = 4351
  Output Stream = Teuchos::RCP<std::__1::basic_ostream<char, std::__1::char_traits<char> > >{ptr=0x7feab0f04ac0,node=0x7feab0d49f20,strong_count=11,weak_count=0}
  Error Stream = Teuchos::RCP<std::__1::basic_ostream<char, std::__1::char_traits<char> > >{ptr=0x7feab0f04ac0,node=0x7feab0d49ff0,strong_count=6,weak_count=0}
 Direction -> 
  Method = Newton
  Newton -> 
   Forcing Term Method = Constant
   Rescue Bad Newton Solve = 1
   Linear Solver -> 
    Max Iterations = 500   [unused]
    Tolerance = 1e-08
    Write Linear System = 0   [unused]
 Linear Solver -> 
  Jacobian Operator = Matrix-Free
  Zero Initial Guess = 1
  Aztec Solver = GMRES
  Output Frequency = 1
  Preconditioner = New Ifpack
  Preconditioner Reuse Policy = Reuse
  Compute Scaling Manually = 0
  Output Solver Details = 1   [default]
  Throw Error on Prec Failure = 1   [default]
  RCM Reordering = Disabled   [default]
  Orthogonalization = Classical   [default]
  Size of Krylov Subspace = 300   [default]
  Convergence Test = r0   [default]
  Max Age Of Prec = 1   [default]
 Line Search -> 
  Method = Backtrack
  Backtrack -> 
   Minimum Step = 1e-12   [default]
   Default Step = 1   [default]
   Recovery Step = 1   [default]
   Max Iters = 100   [default]
   Reduction Factor = 0.5   [default]
 Solver Options -> 
  Status Test Check Type = Minimal   [default]

And parameters used for creating preconditioner for the VERY first nonlinear step:

-- Nonlinear Solver Step 0 -- ||F|| = 3.123e-03 step = 0.000e+00 dx = 0.000e+00

   CALCULATING FORCING TERM
   Method: Constant
   Forcing Term: 1e-08

   Creating a new preconditioner

NOX::Epetra::LinearSolverAztecOO : createNewIfpackPrecon - using Teuchos parameter : [empty list]

   Time required to create preconditioner : 0.0189149 (sec.)

    *******************************************************
    ***** Problem: NOX::Matrix-Free
    ***** Preconditioned GMRES solution
    ***** IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)
    ***** No scaling
    *******************************************************

            iter:    0           residual = 1.000000e+00
            iter:    1           residual = 1.733014e-02

......

After 52 iterations, I got the following warning message:

        Solver:         gmres
number of iterations:   52

Actual residual =  5.4065e-05   Recursive residual =  2.6018e-11

Calculated Norms                Requested Norm
--------------------------------------------    --------------

||r||_2 / ||r0||_2:     1.731280e-02    1.000000e-08

***************************************************************

    Solution time: 0.033199 (sec.)
    total iterations: 52

WARNING: NOX::Direction::Newton::compute() - Linear solve failed to achieve convergence - using the step anyway since "Rescue Bad Newton Solve" is true.

Then the nonlinear step 1 starts, and used the following default Teuchos parameters in creating the preconditioner:

-- Nonlinear Solver Step 1 -- ||F|| = 1.101e-04 step = 1.000e+00 dx = 3.823e-01

   CALCULATING FORCING TERM
   Method: Constant
   Forcing Term: 1e-08

   Preconditioner Reuse: Age of Prec --> 1 / 1

   Destroying preconditioner

   Creating a new preconditioner

NOX::Epetra::LinearSolverAztecOO : createNewIfpackPrecon - using Teuchos parameter : fact: relax value = 0 [default] fact: absolute threshold = 0 [default] fact: relative threshold = 1 [default] fact: level-of-fill = 0 [default]

   Time required to create preconditioner : 0.00630021 (sec.)

    *******************************************************
    ***** Problem: NOX::Matrix-Free
    ***** Preconditioned GMRES solution
    ***** IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)
    ***** No scaling
    *******************************************************
rppawlo commented 6 years ago

As a sanity check, can you fix your preconditioning matrix to identity for the entire run? Then you should get the same result as if preconditioning is off but still run through all the preconditioner machinery.

guoluanjing commented 6 years ago

Will do. But before that, just to make sure the preconditioner machinery has been utilized correctly: do I pass the preconditioning matrix A (an Epetra::CrsMatrix object) into the LinearSystemAztecOO as the EpetraOperator M?

rppawlo commented 6 years ago

Yes - that is correct. At a minimum, it must be an EetraOperator. To use internal preconditioners like ifpack or ml, we need to cast the matrix to an Epetra_RowMatrix. The CrsMatrix derives from that class.

guoluanjing commented 6 years ago

Hi Roger @rppawlo, I did the sanity check. The identity matrix works and outputs results that are essentially the same to the case with no preconditioning. So from here, what do you suggest me do in order to get my actual preconditioning matrix working?

guoluanjing commented 6 years ago

Actually, this is interesting. Is it possible that with an identity preconditioning matrix the numerical performance is better (less nonlinear iterations)?

rppawlo commented 6 years ago

So from here, what do you suggest me do in order to get my actual preconditioning matrix working?

I guess I would need to know some more details about what you are implementing as a preconditioner matrix. As a first step, NOX has a finite different matrix that you can use to generate a Jacobian. Brute force finite differencing is very inefficient, but this can be used on small test problems to verify your matrix implementation. If you are around next week I can come into ORNL to talk some more.

guoluanjing commented 6 years ago

We use property tables to lookup parameter values and derivatives. So the finite difference matrix has a hard time capturing the constitutive relationships represented by the tables.

I would love to meet and talk more next week. My schedule is flexible. Please let me know which day and time would work for you.

Luanjing

Oak Ridge, TN

On Fri, Nov 17, 2017 at 7:56 AM, Roger Pawlowski notifications@github.com wrote:

So from here, what do you suggest me do in order to get my actual preconditioning matrix working?

I guess I would need to know some more details about what you are implementing as a preconditioner matrix. As a first step, NOX has a finite different matrix that you can use to generate a Jacobian. Brute force finite differencing is very inefficient, but this can be used on small test problems to verify your matrix implementation. If you are around next week I can come into ORNL to talk some more.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/trilinos/Trilinos/issues/1991#issuecomment-345236882, or mute the thread https://github.com/notifications/unsubscribe-auth/AXj0jAaqGx0DW9Ni_pNuGzHD4wB0iobMks5s3YKDgaJpZM4QfoAC .