[Trilinos][AMGCL] amgcl_schur_complement not working

[sunethwarna]

Hi everyone,

I ran into an error when using "amgcl_schur_complement" linear solver for fluid dynamics problem. I am using the latest kratos master (7.0.0-91521e42a4-Release). The error I'm getting is: "Error: Error in subdomain_deflation parameters: def_vec is not set"

Attached here with is the error log. Could you please have a look @ddemidov @RiccardoRossi ? I will attach a small test case as well to replicate the error.

`World size: 2
Traceback (most recent call last):
  File "MainKratos.py", line 32, in <module>
    simulation.Run()
  File "/home/suneth/software/kratos_master/kratos/python_scripts/analysis_stage.py", line 50, in Run
    self.RunSolutionLoop()
  File "/home/suneth/software/kratos_master/kratos/python_scripts/analysis_stage.py", line 67, in RunSolutionLoop
    is_converged = self._GetSolver().SolveSolutionStep()
  File "/home/suneth/software/kratos_master/applications/FluidDynamicsApplication/python_scripts/fluid_solver.py", line 103, in SolveSolutionStep
    is_converged = self.solver.SolveSolutionStep()
RuntimeError: Error: Error in subdomain_deflation parameters: def_vec is not set

in /home/suneth/software/kratos_master/applications/TrilinosApplication/external_includes/amgcl_mpi_schur_complement_solver.h:270:bool AmgclMPISchurComplementSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>::Solve(AmgclMPISchurComplementSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>::SparseMatrixType&, AmgclMPISchurComplementSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>::VectorType&, AmgclMPISchurComplementSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>::VectorType&) [with TSparseSpaceType = TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>; TDenseSpaceType = UblasSpace<double, boost::numeric::ublas::matrix<double>, boost::numeric::ublas::vector<double> >; TReordererType = Reorderer<TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>, UblasSpace<double, boost::numeric::ublas::matrix<double>, boost::numeric::ublas::vector<double> > >; AmgclMPISchurComplementSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>::SparseMatrixType = Epetra_FECrsMatrix; AmgclMPISchurComplementSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>::VectorType = Epetra_FEVector]
   /home/suneth/software/kratos_master/applications/TrilinosApplication/custom_strategies/builder_and_solvers/trilinos_block_builder_and_solver.h:365:void TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::SystemSolveWithPhysics(TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemMatrixType&, TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemVectorType&, TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemVectorType&, ModelPart&) [with TSparseSpace = TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>; TDenseSpace = UblasSpace<double, boost::numeric::ublas::matrix<double>, boost::numeric::ublas::vector<double> >; TLinearSolver = LinearSolver<TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>, UblasSpace<double, boost::numeric::ublas::matrix<double>, boost::numeric::ublas::vector<double> > >; TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemMatrixType = Epetra_FECrsMatrix; TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemVectorType = Epetra_FEVector]
   /home/suneth/software/kratos_master/applications/TrilinosApplication/custom_strategies/builder_and_solvers/trilinos_block_builder_and_solver.h:422:void TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::BuildAndSolve(typename TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::BaseType::TSchemeType::Pointer, ModelPart&, TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemMatrixType&, TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemVectorType&, TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemVectorType&) [with TSparseSpace = TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>; TDenseSpace = UblasSpace<double, boost::numeric::ublas::matrix<double>, boost::numeric::ublas::vector<double> >; TLinearSolver = LinearSolver<TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>, UblasSpace<double, boost::numeric::ublas::matrix<double>, boost::numeric::ublas::vector<double> > >; typename TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::BaseType::TSchemeType::Pointer = shared_ptr<Scheme<TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>, UblasSpace<double, boost::numeric::ublas::matrix<double>, boost::numeric::ublas::vector<double> > > >; TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemMatrixType = Epetra_FECrsMatrix; TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::TSystemVectorType = Epetra_FEVector]

Test case files

[ddemidov]

@sunethwarna , it looks like these two lines:

https://github.com/KratosMultiphysics/Kratos/blob/5f38d58d3be7bea856e629c06e8f24d6ec960beb/applications/TrilinosApplication/external_includes/amgcl_mpi_schur_complement_solver.h#L228-L229

should be replaced with

        mprm.put("precond.psolver.num_def_vec", 1);
        mprm.put("precond.psolver.def_vec", &dv);

Also, precond.psolver.solver should be replaced with precond.psolver.isolver here: https://github.com/KratosMultiphysics/Kratos/blob/5f38d58d3be7bea856e629c06e8f24d6ec960beb/applications/TrilinosApplication/external_includes/amgcl_mpi_schur_complement_solver.h#L164-L166

But, @RiccardoRossi , I think it would be even better to replace amgcl::mpi::subdomain_deflation with an mpi amg here. What do you think?

[sunethwarna]

@ddemidov The suggested changes work for me. Thanks :). I will leave it to @ddemidov and @RiccardoRossi to decide on how to implement changes to the master branch.

[RiccardoRossi]

sorry guys, i am out of office, i ll take a look asap.

btw, with the oroposed changes does it work better than normal amg?

[RiccardoRossi]

@ddemidov i am (of course) ok with the changes you propose. I guess that this changed at some point when we updated the amgcl version.

regarding the possibility of using full_amg instead of deflation i am also in principle ok if you believe it will work. My doubt comes as I seem to remember that we were using a matrix-free method in the solution of the pressure complement so i am not 100% sure it will work correctly with full amg. If my comment is wrong then let's just correct this in a PR

[RiccardoRossi]

BTW @sunethwarna is (i think) solving a stiff problem in terms of two turbulence-related variables. Maybe the CPR preconditioner would work good in that context

[ddemidov]

5753 makes the easy job of fixing the parameters as discussed above. I guess we will need to experiment in order to decide if replacing subdomain deflation with amg for the pressure part solver is worth it.