ILU preconditioner and mesh-scale dependency?

[Jaekwang Kim]

Hi 

I had a trouble in solving scalar transport equation using DG elements. 

For example, I was solving for \lambda in the following equation where u is also a numerical solution

The thing is I do not get any problem when my mesh file is scale of O(1) and I have tested my code works fine using the Method of Manufactured solution. 

Yet, when I tried to solve the governing equation at small scale mesh...(I.e. all (x,y) point in the original mesh-file was reduced by (0.01x,0.01y) 

I receive error like this 

An error occurred in line <129> of file </Users/jaekwangkim/Program/dealii-8.5.0/include/deal.II/lac/sparse_ilu.templates.h> in function

    void dealii::SparseILU<double>::initialize(const SparseMatrix<somenumber> &, const dealii::SparseILU::AdditionalData &) [number = double, somenumber = double]

The violated condition was: 

    luval[ia[k]] != 0

Additional information: 

    While computing the ILU decomposition, the algorithm found a zero pivot on the diagonal of row 34. This must stop the ILU algorithm because it means that the matrix for which you try to compute a decomposition is singular.

The problem occurs when the red line of the code is executed...

 template <int dim>

    void StokesProblem<dim>::solve_transport ()

    {

        std::cout << "   -solve transport begins"<< std::endl;

        SolverControl           solver_control (std::pow(10,6), system_rhs.block(2).l2_norm() * pow(10,-4));

        

        unsigned int restart = 500;

        SolverGMRES< Vector<double> >::AdditionalData gmres_additional_data(restart+2);

        SolverGMRES< Vector<double> > solver(solver_control, gmres_additional_data);

      

        //make preconditioner

        SparseILU<double>::AdditionalData additional_data(0,500); // (0 , additional diagonal terms)

        std::cout << "B" << std::endl;

        SparseILU<double> preconditioner;

        

        preconditioner.initialize (system_matrix.block(2,2), additional_data);

        

        std::cout << "A" << std::endl;

        solver.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), preconditioner);

        // constraints.distribute (solution);

    }

    

I wonder now ...

1. why smaller scale mesh only results such error,

2. If possible, is there any choice of better preconditionner for my system other than SparseILU to go around this problem? 

Thanks,

Regards, 

Jaekwang Kim 

[Wolfgang Bangerth]

On 3/13/19 8:18 PM, Jaekwang Kim wrote:
>
>
> I wonder now ...
>
> 1. why smaller scale mesh only results such error,

You have two terms in your matrix, one from the advection term and one
from the reaction term. If you scale the mesh by a factor of c, the
height of the shape functions used in the reaction term does not change,
but the gradient of shape functions used in the advection term is scaled
by 1/c. As a consequence, the balance of the two terms that together
make up the matrix changes and you apparently get into a situation where
the resulting matrix is ill-conditioned or indeed not invertible.

> 2. If possible, is there any choice of better preconditionner for my
> system other than SparseILU to go around this problem?

To make sure that the system you have assembled is correct, I usually
first check with the SparseDirectUMFPACK solver. Once I know that the
system is correct, I think about iterative solvers and preconditioners.
In your case, the system you have is advection dominated, for which it
is notoriously difficult to construct good preconditioners.

Best
WB

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Jaekwang Kim]

Thanks for sharing insight!