Solving an interface problem using matrix-free

[David Schneider]

Hello all together,

I'm currently working on a modified version of step-37 dealing with a matrix-free implementation of a Laplace problem. As a postprocessing step, I need to solve an interface problem: the function is essentially already implemented in VectorTools::project_boundary_values(). However, I cannot use the default implementation in the library, because I don't have a boundary function to project, but I have already a nodal vector (the RHS) describing the boundary.

Thus, I copied more or less the implementation of 'VectorTools::project_boundary_values' into my program and this works out. However, the implementation relies on a matrix-based approach and there are lots of data structures I need to setup for this small and simple system. I was wondering if I can re-use the MatrixFree object for this purpose, as it has already all the preconditioners and operators implemented.

The idea would be to pass a global vector into the cg.solve() function and (instead of looping over all cells) looping over all boundary cells and only perform operations on the desired faces belonging to the relevant interface. There are mainly three questions I was wondering about: (1) will the SolverControl emit convergence although I operate only on some of the DoFs (do I need to reset the remaining ones) ? (2) are there other options to achieve this (solving a problem on a  subdomain such as an interface) using matrix-free (3) do you think the approach makes sense or should I just go with the sparse matrix-vector approach since the overall system will be anyway relatively small (dim - 1)?

Thanks in advance and kind regards,
David

[David Schneider]

I can answer it myself after a bit of work, maybe there is someone wondering about similar things in the future:

> (1) will the SolverControl emit convergence although I operate only on some of the DoFs (do I need to reset the remaining ones)

Yes it will. Just reset the irrelevant DoFs and consider only the relevant here.

> (2) are there other options to achieve this (solving a problem on a  subdomain such as an interface) using matrix-free

Not sure. However, project_boundary_values works only for serial Triangulations (there is actually an Assertion in deal.II missing). Reducing the overall system size to the size of the subdomain (as it is done for the serial version) while still considering the global coupling of distributed DoFs is exceedingly complicated. With matrix-free, I only touch the relevant DoFs when solving the system and I still have the DoF coupling and ghost value exchange as usual. The only 'price' I pay consists of the additional memory consumption for a global vector in the CGSolve, which is tolerable, or at least the best I can get.

> (3) do you think the approach makes sense or should I just go with the sparse matrix-vector approach since the overall system will be anyway relatively small (dim - 1)?

Makes sense. Choosing here sparse matrices is not a good idea as then I would have more unused memory allocated. Also, the (my selected) Jacobi preconditioner might become problematic for sparse matrices, if DoFs on the diagonal are ignored completely.