[dune-pdelab] Fwd: solver fails to reset correctly after FMatrixError (singular matrix)

[Shubhangi Gupta]

Dear pdelab users,

I am currently experiencing a rather strange problem during parallel 
solution of my finite volume code. I have written a short outline of my 
code below for reference.

At some point during computation, if dune throws an error, the code 
catches this error, resets the solution vector to the old value, halves 
the time step size, and tries to redo the calculation (osm.apply()).

However, if I get the error "FMatrixError: matrix is singular", the 
solver seems to freeze. Even the initial defect is not shown! (See the 
terminal output below.) I am not sure why this is so, and I have not 
experienced this issue before.

I will be very thankful if someone can help me figure out a way around 
this problem.

Thanks, and warm wishes, Shubhangi


*// code layout*

     ...UG grid, generated using gmsh, GV, ...

     typedef Dune::PDELab::QkDGLocalFiniteElementMap<GV::Grid::ctype, 
double, 0, dim, Dune::PDELab::QkDGBasisPolynomial::lagrange> FEMP0;
     FEMP0 femp0;
     typedef 
Dune::PDELab::GridFunctionSpace<GV,FEMP0,Dune::PDELab::P0ParallelConstraints,Dune::PDELab::ISTL::VectorBackend<>> 
GFS0;
     GFS0 gfs0(gv,femp0);
     typedef Dune::PDELab::PowerGridFunctionSpace< GFS0,num_of_vars, 
Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>, 
Dune::PDELab::EntityBlockedOrderingTag> GFS_TCH;

     ... LocalOperator LOP lop, TimeLocalOperator TOP top, GridOperator 
GO go, InstationaryGridOperator IGO igo, ...

     typedef Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<IGO> LS;
     LS ls(gfs,50,1,false,true);
     typedef Dune::PDELab::Newton< IGO, LS, U > PDESOLVER;
     PDESOLVER pdesolver( igo, ls );
Dune::PDELab::ImplicitEulerParameter<double> method;

     Dune::PDELab::OneStepMethod< double, IGO, PDESOLVER, U, U > osm( 
method, igo, pdesolver );

     //TIME-LOOP
     while( time < t_END - 1e-8){
             try{
                 //PDE-SOLVE
                 osm.apply( time, dt, uold, unew );
                 exceptionCaught = false;
             }catch ( Dune::Exception &e ) {
                 //RESET
                 exceptionCaught = true;
                 std::cout << "Catched Error, Dune reported error: " << 
e << std::endl;
                 unew = uold;
                 dt *= 0.5;
osm.getPDESolver().discardMatrix();
                 continue;
             }
             uold = unew;
             time += dt;
     }


*// terminal output showing FMatrixError...*


  time = 162.632 , time+dt = 164.603 , opTime = 180 , dt  : 1.97044

  READY FOR NEXT ITERATION.
_____________________________________________________
  current opcount = 2
****************************
TCH HYDRATE:
****************************
TIME STEP [implicit Euler]     89 time (from):   1.6263e+02 dt:   
1.9704e+00 time (to):   1.6460e+02
STAGE 1 time (to):   1.6460e+02.
   Initial defect:   2.1649e-01
Using a direct coarse solver (SuperLU)
Building hierarchy of 2 levels (inclusive coarse solver) took 0.2195 
seconds.
=== BiCGSTABSolver
  12.5        6.599e-11
=== rate=0.1733, T=1.152, TIT=0.09217, IT=12.5
   Newton iteration  1.  New defect:   3.4239e-02.  Reduction (this):   
1.5816e-01.  Reduction (total):   1.5816e-01
Using a direct coarse solver (SuperLU)
Building hierarchy of 2 levels (inclusive coarse solver) took 0.195 seconds.
=== BiCGSTABSolver
    17        2.402e-11
=== rate=0.2894, T=1.655, TIT=0.09738, IT=17
   Newton iteration  2.  New defect:   3.9906e+00.  Reduction (this):   
1.1655e+02.  Reduction (total):   1.8434e+01
Using a direct coarse solver (SuperLU)
Building hierarchy of 2 levels (inclusive coarse solver) took 0.8697 
seconds.
=== BiCGSTABSolver
Catched Error, Dune reported error: FMatrixError 
[luDecomposition:/home/sgupta/dune_2_6/source/dune/dune-common/dune/common/densematrix.hh:909]: 
matrix is singular
_____________________________________________________
  current opcount = 2
****************************
TCH HYDRATE:
****************************
TIME STEP [implicit Euler]     89 time (from):   1.6263e+02 dt:   
9.8522e-01 time (to):   1.6362e+02
STAGE 1 time (to):   1.6362e+02.

*... nothing happens here... the terminal appears to freeze...*



-- 
Dr. Shubhangi Gupta
Marine Geosystems
GEOMAR Helmholtz Center for Ocean Research
Wischhofstraße 1-3,
D-24148 Kiel

Room: 12-206
Phone: +49 431 600-1402
Email:sgupta at geomar.de

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20190710/afe90eb9/attachment.htm>