4
$\begingroup$

I am reading "Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations" (Smith 1996), and I am confused as to how the below Alternating Schwarz algorithm (algorithm 1.1.2 from Smith 1996), specifically the term $v$ would be incorporated into a preconditioned krylov method (e.g., conjugate gradient) to solve the linear system in eq. (1.8). Is $v$ equivalent to the preconditioning matrix $M$? Is it the vector $\vec{z}$ as shown in the conjugate gradient wiki?

Here is the algorithm:

enter image description here

and note that from the book "the application $A_{\Omega_1}^{-1}$ and $A_{\Omega_2}^{-1}$ may also be replaced with some suitable approximate solver". I think this would mean the preconditioning itself would also require something like its own iterative or direct solver (which seems to also be corroborated by p. 148 of ch.6 in "PETSc for Partial Differential Equations" (Bueler 2021)).

Here is the linear system:

enter image description here

Based on the conjugated gradient method preconditioned with the additive Schwarz method described in ch. 3 of "Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation" (Dolean 2015), I would think the $v$ resulting from the alternating Schwarz algorithm shown in my question would have to correspond to $\vec{z}$, and it wouldn't be the matrix $M$ just based on the fact that the matrix-vector multiplications would result in a vector. But I am not 100% sure.

Edit:

In a very abstract sense, I think essentially the following algorithm would be used to solve a given linear system (say arising from a finite element discretization) using a domain decomposition method (in this case Alternating Schwarz) as a preconditioner:

# Set up global linear system Au = f
A, f = discretize(governing_pde) 
ndofs = length(f) # number of degrees of freedom
u = zeros(ndofs)

# Partition domain
domain = partition_domain(A)
internal_domain_omega_indices = get_internal_domain_indices(domain)
boundary_domain_partial_omega_indices = get_boundary_domain_indices(domain)
interface_gamma_indices = get_interface_indices(domain)

# Initialize a solver to be used within the domain decomposition method
direct_solver = ForwardBackwardSubstitution(factorization = :cholesky)

# Krylov method (e.g., preconditioned conjugate gradient)
# to solve global system
for k in niters:
  v = alternating_schwarz_algo(
    A, u, f, 
    internal_domain_omega_indices,
    boundary_domain_partial_omega_indices,
    interface_gamma_indices,
    internal_solver_for_A_inv = direct_solver)

  # Single iteration of conjugate gradient method that
  # updates the solution to the linear system Au = f inplace
  conjugate_gradient_iteration!(
    A = A, f = f, u = u, preconditioner_z = v)
end

. Here is a related question/answer I wrote on stack overflow.

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.