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

[Andrea Mola]

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_EpetraModelEvaluator.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