Stratimikos Package Browser (Single Doxygen Collection): test_stratimikos_belos_nan

Go to the documentation of this file.
00001 
00002 #include "Thyra_BelosLinearOpWithSolveFactory.hpp"
00003 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00004 #include "Thyra_MultiVectorStdOps.hpp"
00005 #include "Thyra_VectorBase.hpp"
00006 #include "Thyra_VectorStdOps.hpp"
00007 #include "Thyra_EpetraLinearOp.hpp"
00008 #include "EpetraExt_readEpetraLinearSystem.h"
00009 #include "Epetra_SerialComm.h"
00010 
00011 #include "Teuchos_UnitTestHarness.hpp"
00012 
00013 namespace Thyra {
00014 
00015   /* This test is to make sure that stratimikos wrappers to belos
00016      correctly catch belos exceptions thrown and report back a failed
00017      solve.  If belos detects a nan in the linear solve, it will throw
00018      an exception.  This was causing the termination of the
00019      simulation.  We have codes where cutting the time step allows
00020      recovery so it makes sense that stratimikos should catch the
00021      thrown exception and report that the solve failed.  Then the time
00022      integrator can retake a time step.  This test makes sure the
00023      stratimikos interface catches the exception and reports a failed
00024      solve.
00025   */
00026   TEUCHOS_UNIT_TEST(belos, nan_handling)
00027   {
00028     
00029     Epetra_SerialComm comm;
00030     Teuchos::RCP<Epetra_CrsMatrix> epetra_A;
00031     std::string matrixFile = "FourByFour.mtx";
00032     EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
00033 
00034     Teuchos::RCP<const LinearOpBase<double> > A = epetraLinearOp(epetra_A);
00035 
00036     Teuchos::RCP<LinearOpWithSolveFactoryBase<double> >
00037       lowsFactory;
00038     {
00039       Teuchos::RCP<BelosLinearOpWithSolveFactory<double> >
00040         belosLowsFactory = Teuchos::rcp(new BelosLinearOpWithSolveFactory<double>());
00041       lowsFactory = belosLowsFactory;
00042     }
00043 
00044     Teuchos::ParameterList belosLOWSFPL;
00045     {
00046       belosLOWSFPL.set("Solver Type","Block GMRES");
00047 
00048       Teuchos::ParameterList& belosLOWSFPL_solver =
00049   belosLOWSFPL.sublist("Solver Types");
00050       
00051       Teuchos::ParameterList& belosLOWSFPL_gmres =
00052   belosLOWSFPL_solver.sublist("Block GMRES");
00053       
00054       belosLOWSFPL_gmres.set("Maximum Iterations",int(4));
00055       belosLOWSFPL_gmres.set("Convergence Tolerance",double(1.0e-4));
00056       belosLOWSFPL_gmres.set("Maximum Restarts",int(0));
00057       belosLOWSFPL_gmres.set("Block Size",int(1));
00058       belosLOWSFPL_gmres.set("Num Blocks",int(4));
00059       belosLOWSFPL_gmres.set("Output Frequency",int(1));
00060       belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(false));
00061       
00062       lowsFactory->setParameterList(Teuchos::rcp(&belosLOWSFPL,false));
00063     }
00064 
00065     Teuchos::RCP<LinearOpWithSolveBase<double> > nsA = lowsFactory->createOp();
00066     Thyra::initializeOp<double>(*lowsFactory,  A, nsA.ptr());
00067 
00068     Teuchos::RCP<Thyra::VectorBase<double> > x = Thyra::createMember(A->domain());
00069     Teuchos::RCP<Thyra::VectorBase<double> > f = Thyra::createMember(A->range());
00070 
00071     Thyra::put_scalar(0.0,x.ptr());
00072     Thyra::put_scalar(1.0,f.ptr());
00073     
00074     // Insert a nan
00075     Thyra::set_ele(0,Teuchos::ScalarTraits<double>::nan(),x.ptr());
00076 
00077     Thyra::SolveStatus<double> status = Thyra::solve<double>(*nsA,Thyra::NOTRANS,*f,x.ptr());
00078 
00079     TEST_EQUALITY(status.solveStatus, Thyra::SOLVE_STATUS_UNCONVERGED);
00080   }
00081   
00082 }