[aspect-devel] Direct Solver in Aspect

[Timo Heister]

>> No, I plan to use TrilinosWrappers::SolverDirect that can only deal
>> with a Trilinos SparseMatrix (not a Trilinos::BlockSparseMatrix).
>
> How difficult would it be to extend TrilinosWrappers::SolverDirect?

It can not possibly work with a 2x2 table of Trilinos matrices. The
only option I see would be to copy everything from a block matrix into
a single matrix. I don't think this is a good solution.

>>>> 2. fix all the places that assume the current layout (say T is in block
>>>> 3)
>>>
>>> Can this be done through appropriate modifications/additions to the
>>> Introspection class, rather than scattering if/elses throughout?
>>
>> As much as possible yes. But there are many places where we assume
>> certain block indices (I am going through them right now).
>
> It would be a shame if we needed to push this kind of input flag through all
> sorts of code. I hope it can be done without too much clutter in places that
> really don't need to know this :-(

I only need an if in very few places so far. It mostly involves
replacing things like

  IndexSet range = this->get_solution().block(2).locally_owned_elements();

by

  const unsigned int blockidx = this->introspection().block_indices.temperature;
  IndexSet range =
this->get_solution().block(blockidx).locally_owned_elements();

which is useful anyway.

> Anyway, good feature. I'm looking forward to seeing how it works in practice
> :-)

One of the things that I got out of the CIG MC meeting that might help
with the speed for large 2d computations on a single workstation.
We'll see.

-- 
Timo Heister
http://www.math.clemson.edu/~heister/