[Trilinos-Users] [EXTERNAL] Running Rythmos DAE with a user defined Epetra_Operator as Jacobian

Thu May 22 13:23:24 MDT 2014

Andrea,

If you define the OutArg "W" (which is the full linear solver) and not the OutArgs "W_op" or "W_prec", then that should ensure that a default Stratimikos linear solver and preconditioner will not get used.  It will either then just magically work or it will break (and throw pretty fast I would guess).  Otherwise, I would have to see exactly how you are setting all of this up.  At some level, the Rythmos algorithmic classes have to use the Thyra::LinearOpWithSolveBase object you define as "W".  I suspect that if you define "W_op" then there may be some logic that ignores "W" and/or "W_prec".

-Ross

> -----Original Message-----
> From: trilinos-users-bounces at software.sandia.gov [mailto:trilinos-users-
> bounces at software.sandia.gov] On Behalf Of Andrea Mola
> Sent: Thursday, May 22, 2014 2:36 PM
> To: trilinos-users at software.sandia.gov
> Subject: Re: [Trilinos-Users] [EXTERNAL] Running Rythmos DAE with a user
> defined Epetra_Operator as Jacobian
> 
> Dear all,
> I am continuing to struggle in implementing a Rythmos DAE. As previously
> suggested in this thread, I have converted the EpetraExt::ModelEvaluator
> to a Thyra::ModelEvaluator, and managed to run Rythmos with a
> user-defined Epetra_Operator as Jacobian.
> Unfortunately, I still haven't had success in using a used-defined
> Epetra_Operator as Jacobian preconditioner.
> Within my Thyra::ModelEvaluator, I define a Thyra::linearOpWithSolve and
> initialize it with the initializePreconditionedOp method so as to
> specify my own preconditioner. If I try a "test run" for the resulting
> linear system with solve, it works correctly and the preconditioner
> matrix-vector (the Apply member of the preconditioner Epetra_Operator)
> is correctly called multiple times during the system resolution.
> But for some reason, during the time integration Rythmos is completely
> ignoring the preconditioner, and the corresponding matrix-vector product
> is never called. The OutArgs of my Thyra::ModelEvaluator support W,
> W_op, and W_prec. I thought that maybe I should re-initialize
> linearOpWithSolve at each evaluation. But this is not even possible,
> since the W and W_prec pointers assigned during the time integration to
> the evalModel function are always null.
> I initially thought that Eric's concerns regarding the user defined
> preconditioner in Rythmos were only referred to
> EpetraExt::ModelEvaluator, but I am now afraid that they were also
> referred to the Thyra class. So is there a way to provide a user-defined
> preconditioner to Rythmos? Has anyone succeeded in doing that?
> Thanks a lot for the help.
> 
> 
>             Andrea
> 
> 
> 
> 
> On 2014-05-05 18:34, Andrea Mola wrote:
> > Eric,
> > I am currently using version 11.2.3, but maybe the issue isn't there.
> > I realized the problem was most likely my fault. The create_W method
> > in the model evaluator class was returning an Epetra_CrsMatrix. I
> > never called that method, so I din't know it was effectively used.
> > Nonetheless, I now run into a further exception
> >
> > terminate called after throwing an instance of 'std::logic_error'
> >   what():
> > /home/amola/trilinos-11.2.3-
> Source/packages/thyra/adapters/epetraext/src/model_evaluator/Thyra_Ep
> etraModelEvaluator.cpp:751:
> >
> > Throw number = 1
> >
> > Throw test that evaluated to true: !allFuncsWhereScaled
> >
> > Error, we can not currently handle epetra output objects that could
> > not be scaled.  Special code will have to be added to handle this
> > (i.e. using implicit diagonal and multiplied linear operators to
> > implicitly do the scaling.
> >
> > which is maybe called because in my operator some methods are not
> > properly set.
> > However, given the possible preconditioning problems you are hinting
> > at, I'm converting everything into a Thyra model evaluator, as
> > previously suggested by Roscoe.
> > Thanks for the help.
> >
> >
> >        Andrea
> >
> >
> >
> > On 2014-05-05 16:36, Phipps, Eric T wrote:
> >> Andrea,
> >>
> >> What version of Trilinos are you using?  I believe Andy Salinger
> >> recently
> >> fixed some of these issues when using an Epetra_Operator generated
> >> from an
> >> EpetraExt::ModelEvaluator in Rythmos.  However I believe those fixes
> >> are
> >> only in the development branch which should be accessible from the
> >> public
> >> git repository.  Even with the fix for the failed dynamic casts,
> >> there may still be an issue with supplying a user-defined
> >> preconditioner
> >> to Rythmos.
> >>
> >> -Eric
> >>
> >> On 5/5/14, 7:45 AM, "Andrea Mola" <amola at sissa.it> wrote:
> >>
> >>> Dear all,
> >>> I am trying to use Ryrhmos to solve a DAE system for the time
> >>> advancement of a PDE system.
> >>> So far, I have already been able to set up a small example which
> >>> solves
> >>> correctly a small DAE using a Jacobian defined as an
> >>> Epetra_CrsMatrix.
> >>> Unfortunately, I need to use a different Linear Operator, which I
> >>> define
> >>> as an Epetra_Operator (JacobianOperator :public Epetra_Operator). The
> >>> problem comes when I define my AztecOOLinearOpWithSolveFactory to
> be
> >>> passed to my EpetraExt::ModelEvaluator constructor. Even though I try
> >>> to
> >>> specify through the Thyra::initializeOp function that I will use a
> >>> JacobianOperator instead of an Epetra_CrsMatrix, it still seems to
> >>> think
> >>> that the jacobian of my system is an Epetra_CrsMatrix. In fact, upon
> >>> execution, when the model evaluator calls for a Jacobian evaluation,
> >>> the
> >>> dynamic cast
> >>>
> >>> JacobianOperator &jac = dyn_cast<JacobianOperator>(*W_out);
> >>>
> >>> fails with the following exception
> >>>
> >>> terminate called after throwing an instance of 'Teuchos::m_bad_cast'
> >>>   what():
> >>> /home/amola/trilinos-11.2.3-
> Source/packages/teuchos/core/src/Teuchos_dyn_c
> >>> ast.cpp:61:
> >>>
> >>> Throw number = 1
> >>>
> >>> Throw test that evaluated to true: true
> >>>
> >>> dyn_cast<16JacobianOperator>(15Epetra_Operator) : Error, the object
> >>> with
> >>> the concrete type '16Epetra_CrsMatrix' (passed in through the
> >>> interface
> >>> type '15Epetra_Operator')  does not support the interface
> >>> '16JacobianOperator' and the dynamic cast failed!
> >>> Aborted.
> >>>
> >>> So, my question is: what should I do to be able to use an arbitrary
> >>> Epetra_Operator in an EpetraExt::ModelEvaluator? I am not very
> >>> familiar
> >>> with the parameter lists, so probably the answer is there. Or,
> >>> Thyra::initializeOp is just the wrong function to do what I need. Can
> >>> anyone point me out an example that does that, or suggest me some
> >>> solution?
> >>> Thank you very much in advance for the help.
> >>> Regards,
> >>>
> >>>
> >>>     Andrea
> >>> _______________________________________________
> >>> Trilinos-Users mailing list
> >>> Trilinos-Users at software.sandia.gov
> >>> https://software.sandia.gov/mailman/listinfo/trilinos-users
> > _______________________________________________
> > Trilinos-Users mailing list
> > Trilinos-Users at software.sandia.gov
> > https://software.sandia.gov/mailman/listinfo/trilinos-users
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at software.sandia.gov
> https://software.sandia.gov/mailman/listinfo/trilinos-users