fdrmrc / Polydeal

C++ implementation of Polygonal Discontinuous Galerkin method within the deal.II Finite Element library.
https://fdrmrc.github.io/Polydeal/
Other
0 stars 0 forks source link

Transfer operator for agglomeration #98

Closed fdrmrc closed 3 months ago

fdrmrc commented 8 months ago

Depends on #97, first relevant commit is 3fe73d0dff89903e97a4310631ccf81e006b246b

This PR adds a class MGAgglomerationTransfer, which derives form dealii::MGTransferBase. Its reinit() creates the transfer operator (prolongator) from a coarser space to a finer one. Usual prolongate() and prolongate_and_add() simply performs mat-vec products.

Tested so far with a structured square grid and an external mesh (4367a765b4cb13ed43a76cce4cf5ccc2694e49d2), where levels are generated through the rtree. The test interpolates onto the coarse DG space a known (constant and linear) function $f$, call it $u_c$, and prolongates it to the finer grid. After that, the difference between the prolonged vector $u_f$ and the interpolant of $f$ onto the finer space is computed. In the image below, points are the nodal values of the prolonged vector $u_f$, lying on the surface of $u_c$.

Here's the result of the test

N elements (coarse) = 4
N elements (fine) = 16
Construct transfer operator
f= 1
Norm of error(L2): 0
f= x+y
Norm of error(L2): 0
f= x
Norm of error(L2): 0
f= y
Norm of error(L2): 0

immagine

fdrmrc commented 8 months ago

@luca-heltai With this in place I should be able to setup a multigrid test. I'm working on it.

luca-heltai commented 8 months ago

This is very specific to agglomeration handler, in the way you wrote it. Now that you put your hands into it, wouldn't it be better to have a class MGLinearOperatorTransfer, that uses linear operators for the transfers (i.e., the vmult and vmult_add of a linear operator), and one free function that creates a LinearOperator that performs the transfer? This separates the concerns. The class you created needs to know too much about the underlying problem. The transfer is really just a LinearOperator, so there should be some way to build one by passing the two const AgglomerationHandler<dim> &agglo_handler_fine, const AgglomerationHandler<dim> &agglo_handler_coarse, i.e.,

LinearOpearator<VectorType> transfer = transfer_operator(agglo_fine, agglo_coarse);
MGLinearOperatorTransfer<VectorType> mg_transfer(transfer);

This would make MGLinearOperatorTransfer usable for all our other purposes.

luca-heltai commented 8 months ago

Btw, the same is true for #97: build a MGLinearOperatorCoarseSolver such as

Ainv = inverse_operator(A, solver, preconditioner);
MGLinearOperatorCoarseSolver<VectorType> mg_coarse(Ainv);

leave it to the outside world to build the LinearOperator, and have the MG* classes simply interpret externally provided LinearOperators.

fdrmrc commented 8 months ago

CI is failing due to a private member I was accessing while debugging. The latest version has to be cleaned a lot, but I pushed it here for the time being.

fdrmrc commented 3 months ago

Superseded by #123