How to use external linear solvers and preconditioners

From KratosWiki
Jump to: navigation, search

The Kratos shared-memory version, implements an in-house implementation of some iterative solvers as well as interfaces to various external linear solvers and preconditioners. The rationale at the basis of such choice, is that while iterative solvers can be readily implemented, preconditioners and direct solvers are much more technical and difficult. Hence the choice was made to make available state-of-the art libraries developed outside of the kratos.

Such external libraries are made available as either part of the ExternalSolversApplication or of the MKLSolversApplication which need to be included to allow the use of the external solvers.

Additionally, and for the specific case of OpenCL capable sistems, other solvers are available as a part of the OpenCLApplication

Contents

Solvers in the core Kratos

The core of the Kratos natively implements some of the most common Iterative Solvers available in the literature, as well as some preconditioners. Such solvers are directly available once the kratos core is loaded.

CG solver

This is a parallel implementation of the standard CG solver. All of the operations performed within the solver are OpenMP parallel. Solver is suitable for symmetric positive-definite matrices.

The solver can be constructed as (note that precond is an optional argument with the preconditioner to be used)

   tol = 1e-6 #stopping tolerance
   max_it = 500 #max number of  iterations
   solver = CGSolver(tol,max_it, precond)

BiCGStab solver

This solver implements the classical BiCGStab solver, suitable for non-symmetric matrices. All of the operations performed within the solver are OpenMP parallel.

The solver can be constructed as (note that precond is an optional argument with the preconditioner to be used)

   tol = 1e-6 #stopping tolerance
   max_it = 500 #max number of  iterations
   solver = BICGSTABSolver(tol,max_it, precond)

Deflated CG solver

This solver combined the classical CG solver with a subdomain deflation technique. The solver is OpenMP parallel. Solver is suitable for symmetric positive-definite matrices, it is useful for elongated domains since it implements an approach which is somewhat similar to an AMG technique

The solver can be constructed as

   assume_constant_structure = True #set it to false if the matrix shape changes at every iteration (setting it to true increases the cost if the matrix structure is fixed)
   max_reduced_size = 1000 #max size of the coarse system to be solved internally
   tol = 1e-6 #stopping tolerance
   max_it = 500 #max number of  iterations
   solver = DeflatedCGSolver(tol,max_it,assume_constant_structure,max_reduced_size)

MixedUPSolver

This solver implements an block-preconditioners for the non-symmetic velocity-pressure systems stemming out of the discretization of mixed problems (typically in CFD) It combines a linear solver for the velocity block and one for the pressure block. Such new solver need to be predefined as normal linear solvers

Note that the pressure sub-block will be symmetric only if the A12 and A21 blocks of the original matrix are symmetric

  velocity_linear_solver = velocity_linear_solver
  pressure_linear_solver = pressure_linear_solver
  m = 10 #GMRES krylov space dimension
  tol = 1e-6
  max_it = 10
  linear_solver = MixedUPLinearSolver(velocity_linear_solver,pressure_linear_solver,tol,max_it,m)

an important feature of this solver is that it can be used also by setting

  max_it = 1

in such case any solve becomes similar to one predictor-corrector iteration

Solvers in ExternalSolversApplication

Using such solvers requires the users to import the ExternalSolversApplication by writing in the main python file:

  from KratosMultiphysics.ExternalSolversApplication import *

Once this is done the following solvers are available:

SuperLU Direct Solver

This solver is a freely available sparse direct solver, available at the address [1]. The same address also contains a detailed user manual. It can be used to solve a general sparse non-symmetric matrix. The solver is NOT smp parallel.

In order to define a solver object of this type one should write

   solver = SuperLUSolver()

SuperLU Iterative Solver

The SuperLU solver described just above also comes in preconditioner version as a part of the SuperLU library. It implements an efficient supernodal ILUTP algorithm, suitable for general non-symmetric matrices. The kratos packs it together with an implementation of a gmres solver so that its usage is a little simpler.

If default settings are accepted one can use it by writing

   solver = SuperLUIterativeSolver()

Alternatively the user is allowed to set basic settings using the extended interface. The syntax in this case is as follows:

   tol = 1e-6 #stopping tolerance
   max_it = 300 #max number of GMRES iterations
   restart=150 #size of the GMRES basis
   DropTol=1e-4 #dropping parameters ... see manual 
   FillTol = 1e-2 #dropping parameters ... see manual 
   FillFactor=30 #maximum storage
   solver  =  SuperLUIterativeSolver(tol,max_it,restart,DropTol,FillTol,FillFactor)

ARMS Solver (Algebraic Recursive Multilevel Solver)

This is an implementation of the ARMS solver by Yousef Saad, which comes as a part of the ITSOL library. In kratos it works together with a GMRES solver and works for general non-symmetric matrices

Default settings are availabe by constructing the solver as

   solver = ITSOL_ARMS_Solver()

some user setting can be prescribed by using the constructor as follows:

   tol = 1e-6 #stopping tolerance
   max_it = 300 #max number of GMRES iterations
   restart=150 #size of the GMRES basis
   solver = ITSOL_ARMS_Solver(tol,max_it,restart)

AMGCL Solver ( Algebraic Multi Grid )

The library developed by Denis Demidov , available at the address [2], is a powerfull solver with excellent preconditioners.

   smoother = AMGCLSmoother.ILU0  # or DAMPED_JACOBI or SPAI0 or GAUSS_SEIDEL 
   solver_type = AMGCLIterativeSolverType.CG # or BICGSTAB or GMRES or BICGSTAB_WITH_GMRES_FALLBACK
   verbosity = 1
   max_iterations = 1000
   gmres_size = 50
   tol = 1e-5
   solver =  AMGCLSolver(smoother,solver_type,tol,max_iterations,verbosity,gmres_size)

Solvers in OpenCLApplication

If you have a relatevely powerful graphics card ( a GOOD gaming card, i.e. GTX580), then it might be worth a try to use the GPU to solve the system. Despite the system will be solved using OpenCl code, you will not have to modify anything in your code since the assembly of the system will be done in the standard way, using the CPU. Once the complete system of equations is ready to be solved, it will be passed to the GPU and solved there.

The fastest strategy is using the AMG preconditioners, implemented by Denis Demidov. For large domains the GPU can solve the system in half the time compared to the CPU OpenMP implementation. Have in mind that for small system the transfer time (CPU->GPU) might take longer than the solution itself. Also, large systems might not fit into the GPU Ram (you need approx 1GB Ram for every 1 million elemenents). The OpenMP implementation of the AMG solvers is also really fast, so it is a good idea to test it too since the fastest is case-dependent.

Using such solvers requires the users to import the OpenCLApplication by writing in the main python file:

  from KratosMultiphysics.OpenCLApplication import *

Once this is done the solver can be initialized as:

    tol=1e-5
    max_iterations=1000
    #always keep precision in double, single is a bad idea
    solver_type = OpenCLSolverType.CG # or BiCGStab or GMRES
    preconditioner = OpenCLPreconditionerType.AMG_DAMPED_JACOBI # or AMG_SPAI0 or ILU or NoPreconditioner 
    solver= ViennaCLSolver(tol,10000,OpenCLPrecision.Double,OpenCLSolverType.CG,OpenCLPreconditionerType.AMG_DAMPED_JACOBI)

solvers which come as a part of the MKLSolversApplication

The Kratos has an interface to the (closed) MKL solver. User should have a license of the MKL in order to be able to use this solver. Using the solvers require the user to import the MKLSolversApplication

   from KratosMultiphysics.MKLSolversApplication import *

The interface to the Parallel Direct solver PARDISO, which is a parallel sparse solver for general non-symmetric matrices can be constructed as

   solver = MKLPardisoSolver()
Personal tools
Categories