[petsc-users] Jacobian matrix for dual porosity model

[Barry Smith]

> On Jun 13, 2017, at 10:06 AM, Jed Brown <jed at jedbrown.org> wrote:
> 
> Adrian Croucher <a.croucher at auckland.ac.nz> writes:
> 
>> One way might be to form the whole Jacobian but somehow use a modified 
>> KSP solve which would implement the reduction process, do a KSP solve on 
>> the reduced system of size n, and finally back-substitute to find the 
>> unknowns in the matrix rock cells.
> 
> You can do this with PCFieldSplit type Schur, but it's a lot heavier
> than you might like.

   Is it clear that it would produce much overhead compared to doing a custom "reduction to a smaller problem". Perhaps he should start with this and then profiling can show if there are any likely benefits to "specializing more"?

  Barry

> 
>> Another way might be to form only the reduced-size Jacobian and the 
>> other block-diagonal matrices separately, use KSP to solve the reduced 
>> system but first incorporate the reduction process into the Jacobian 
>> calculation routine, and somewhere a post-solve step to back-substitute 
>> for the unknowns in the matrix cells. However currently we are using 
>> finite differences to compute these Jacobians and it seems to me it 
>> would be messy to try to do that separately for each of the 
>> sub-matrices. Doing it the first way above would avoid all that.
> 
> If you choose this option, you would make your residual evaluation
> perform the local solve (i.e., eliminating the local variables).
> 
>> Any suggestions for what might be a good approach? or any other ideas 
>> that could be easier to implement with PETSc but have similar 
>> efficiency? I didn't see anything currently in PETSc specifically for 
>> solving block-tridiagonal systems.
> 
> Any incomplete or complete factorization (optionally inside block
> Jacobi) is an O(N) direct solver for a tridiagonal system.