Add LST preconditioner

[jowezarek]

We implement preconditioners for mass matrices proposed by Loli, Sangalli & Tani in their paper Easy and efficient preconditioning of the isogeometric mass matrix. Fixes #525.

Obtaining LST preconditioner

Example 1: Working with a DiscreteDeRham object.

# in 3D
M0_pc, M1_pc, M2_pc, M3_pc = derham_h.LST_preconditioner(M0=M0, M1=M1, M2=M2, M3=M3)

returns preconditioner for all four mass matrices. Obtaining only a subset works as follows:

# in 3D
M0_pc, M2_pc = derham_h.LST_preconditioner(M0=M0, M2=M2)

When working with modified mass matrices

# in 3D
DP0, DP1, DP2, _ = derham_h.dirichlet_projectors(kind='linop')
I0, I1, I2       = [IdentityOperator(Vh.coeff_space) for Vh in derham_h.spaces[:-1]]

M0_0 = DP0 @ M0 @ DP0 + (I0 - DP0)
M1_0 = DP1 @ M1 @ DP1 + (I1 - DP1)
M2_0 = DP2 @ M2 @ DP2 + (I2 - DP2)

set the hom_bc flag to True in order to obtain the correct pcs:

# M3_pc remains unchanged for obvious reasons
M0_0_pc, M1_0_pc, M2_0_pc, M3_pc = derham_h.LST_preconditioners(M0=M0, M1=M1, M2=M2, M3=M3, hom_bc=True)

Use these preconditioner like so:

M0_inv   = inverse(M0,   'pcg', pc=M0_pc)
M0_0_inv = inverse(M0_0, 'pcg', pc=M0_0_pc)
M1_0_inv = inverse(M1_0, 'pcg', pc=M1_0_pc)
M2_0_inv = inverse(M2_0, 'pcg', pc=M2_0_pc)
M3_inv   = inverse(M3,   'pcg', pc=M3_pc)

Example 2: Working with FEM objects only. Use the following function located in linalg/tests/test_solvers.py.
It returns a single LST preconditioner for the mass matrix provided.

construct_LST_preconditioner(M, domain_h, fem_space, hom_bc=False)

Use like so

# Given a mass matrix M, a discrete domain domain_h, and a fem_space Vh (corresponding to M)
M_pc    = construct_LST_preconditioner(M, domain_h, Vh)
M_inv   = inverse(M, 'pcg', pc=M_pc, tol=tol, maxiter=maxiter)

DP      = DirichletProjector(Vh)
I       = IdentityOperator(Vh.coeff_space)
M_0     = DP @ M @ DP + (I - DP)
M_0_pc  = construct_LST_preconditioner(M, domain_h, Vh, hom_bc=True)
M_0_inv = inverse(M_0, 'pcg', pc=M_0_pc, tol=tol, maxiter=maxiter)

For more information, check out the new example in examples/vector_potential_3d.py.

Also, check out the test linalg/tests/test_solvers.py::test_LST_preconditioner for more information regarding their efficiency.

[codacy-production]

Coverage summary from Codacy

See diff coverage on Codacy

Coverage variation	Diff coverage
✅ +0.25%	✅ 92.26%

Coverage variation details

	Coverable lines	Covered lines	Coverage
Common ancestor commit (a89a6ee)	31852	19688	61.81%
Head commit (98236b2)	64034 (+32182)	39740 (+20052)	62.06% (+0.25%)

Coverage variation is the difference between the coverage for the head and common ancestor commits of the pull request branch: <coverage of head commit> - <coverage of common ancestor commit>

Diff coverage details

	Coverable lines	Covered lines	Diff coverage
Pull request (#538)	168	155	92.26%

Diff coverage is the percentage of lines that are covered by tests out of the coverable lines that the pull request added or modified: <covered lines added or modified>/<coverable lines added or modified> * 100%

See your quality gate settings    Change summary preferences

[yguclu]

Good job @jowezarek!

I am posting a few comments below

[yguclu]

Good job! I have just some minor comments.

Could you please add the license header to every new file?

[yguclu]

Looks good to me now!