ZeroDimensionalRythmosSolverInterface.cpp 5.39 KB
Newer Older
1
2
//@HEADER
/*
Cristiano Malossi's avatar
Cristiano Malossi committed
3
*******************************************************************************
4

Cristiano Malossi's avatar
Cristiano Malossi committed
5
6
    Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
    Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7

Cristiano Malossi's avatar
Cristiano Malossi committed
8
    This file is part of LifeV.
9

Cristiano Malossi's avatar
Cristiano Malossi committed
10
11
12
13
    LifeV is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
14

Cristiano Malossi's avatar
Cristiano Malossi committed
15
16
17
18
    LifeV is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.
19

Cristiano Malossi's avatar
Cristiano Malossi committed
20
21
    You should have received a copy of the GNU Lesser General Public License
    along with LifeV.  If not, see <http://www.gnu.org/licenses/>.
22

Cristiano Malossi's avatar
Cristiano Malossi committed
23
24
*******************************************************************************
*/
25
26
27
28
//@HEADER

/*!
 *  @file
Cristiano Malossi's avatar
Cristiano Malossi committed
29
 *  @brief Rythmos solver Interface.
30
31
32
33
 *
 *  @date 21-11-2011
 *  @author Mahmoud Jafargholi
 *
34
35
 *  @contributors Cristiano Malossi <cristiano.malossi@epfl.ch>
 *  @mantainer    Cristiano Malossi <cristiano.malossi@epfl.ch>
36
37
 */

38
#include <lifev/zero_dimensional/solver/ZeroDimensionalRythmosSolverInterface.hpp>
39

40
namespace LifeV {
41

42
#if ( defined(HAVE_NOX_THYRA) && defined(HAVE_TRILINOS_RYTHMOS) )
43
44
45
// ===================================================
// Constructors
// ===================================================
46
RythmosSolverInterface::RythmosSolverInterface( Int numCircuitElements,
47
48
                                                Teuchos::RCP< Epetra_Comm > &epetra_comm_ptr,
                                                rythmosModelInterfacePtrRCP_Type theModel ) :
49
    M_epetraCommPtr( epetra_comm_ptr ), M_numElements( numCircuitElements ), M_problemInterfacePtr( theModel ), M_comm( epetra_comm_ptr )
50
{
51
    initialize();
52
53
}

54
void RythmosSolverInterface::initialize()
55
{
56
57
    M_epetraMapPtr = Teuchos::rcp( new Epetra_Map( M_problemInterfacePtr->getMap() ) );
    M_Wgraph = Teuchos::rcp( new Epetra_CrsGraph( M_problemInterfacePtr->getGraph() ) );
58
59
60
}

// Overridden from EpetraExt::ModelEvaluator
61
Teuchos::RCP< const Epetra_Map > RythmosSolverInterface::get_x_map() const
62
{
63
    return M_epetraMapPtr;
64
65
}

66
Teuchos::RCP< const Epetra_Map > RythmosSolverInterface::get_f_map() const
67
{
68
    return M_epetraMapPtr;
69
70
}

71
Teuchos::RCP< const Epetra_Vector > RythmosSolverInterface::get_x_init() const
72
73
{
    Epetra_Vector& soln = M_problemInterfacePtr->getSolutionY();
74
    Teuchos::RCP< Epetra_Vector > x_init = Teuchos::rcp( new Epetra_Vector( soln ) );
75
76
77
    return x_init;
}

78
Teuchos::RCP< const Epetra_Vector > RythmosSolverInterface::get_x_dot_init() const
79
{
80
81
82
    Epetra_Vector& soln = M_problemInterfacePtr->getSolutionYp();
    Teuchos::RCP< Epetra_Vector > x_dot_init = Teuchos::rcp( new Epetra_Vector( soln ) );
    return x_dot_init;
83
84
}

85
Teuchos::RCP< Epetra_Operator > RythmosSolverInterface::create_W() const
86
{
87
88
89
    Teuchos::RCP< Epetra_Operator > W = Teuchos::rcp( new Epetra_CrsMatrix( ::Copy,
                                                                            *M_Wgraph ) );
    return W;
90
91
}

92
EpetraExt::ModelEvaluator::InArgs RythmosSolverInterface::createInArgs() const
93
{
94
    InArgsSetup inArgs;
95
96
97
98
99
    inArgs.setSupports( IN_ARG_x, true );
    inArgs.setSupports( IN_ARG_x_dot, true );
    inArgs.setSupports( IN_ARG_alpha, true );
    inArgs.setSupports( IN_ARG_beta, true );
    inArgs.setSupports( IN_ARG_t, true );
100
    return inArgs;
101
102
}

103
EpetraExt::ModelEvaluator::OutArgs RythmosSolverInterface::createOutArgs() const
104
{
105
    OutArgsSetup outArgs;
106
107
108
109
    outArgs.setSupports( OUT_ARG_f, true );
    outArgs.setSupports( OUT_ARG_W, true );
    outArgs.setSupports( OUT_ARG_W, true );
    outArgs.set_W_properties( DerivativeProperties( DERIV_LINEARITY_NONCONST, DERIV_RANK_UNKNOWN, true ) );
110
    return outArgs;
111
112
}

113
114
void RythmosSolverInterface::evalModel( const InArgs& inArgs,
                                        const OutArgs& outArgs ) const
115
{
116
117
118
119
    Teuchos::RCP< const Epetra_Vector > x = inArgs.get_x();
    Teuchos::RCP< const Epetra_Vector > xdot = inArgs.get_x_dot();
    Real t = inArgs.get_t();

120
#ifdef HAVE_LIFEV_DEBUG
121
122
123
124
125
    std::cout << "RythmosSolverInterface::evalModel ---------------------------{" << std::endl;
    std::cout << "x = " << std::endl;
    x->Print(std::cout);
    std::cout << "xdot = " << std::endl;
    xdot->Print(std::cout);
126
#endif // HAVE_LIFEV_DEBUG
127
128
129
130

    Teuchos::RCP< Epetra_Vector > f;
    if ( ( f = outArgs.get_f() ).get() )
    {
131
        M_problemInterfacePtr->evaluateFImplicit( t, &*x, &*xdot, &*f );
132
#ifdef HAVE_LIFEV_DEBUG
133
134
        std::cout << "f = " << std::endl;
        f->Print(std::cout);
135
#endif // HAVE_LIFEV_DEBUG
136
137
138
139
140
141
142
143

    }
    Teuchos::RCP< Epetra_Operator > W;
    if ( ( W = outArgs.get_W() ).get() )
    {
        const Real alpha = inArgs.get_alpha();
        const Real beta = inArgs.get_beta();
        Epetra_CrsMatrix& jac = Teuchos::dyn_cast< Epetra_CrsMatrix >( *W );
144
        M_problemInterfacePtr->evaluateWImplicit( t, alpha, beta, &*x, &*xdot, &jac );
145
#ifdef HAVE_LIFEV_DEBUG
146
147
        std::cout << "jac = " << std::endl;
        jac.Print(std::cout);
148
#endif // HAVE_LIFEV_DEBUG
149
    }
150
151
152
153
154
#ifdef HAVE_LIFEV_DEBUG
    std::cout << "RythmosSolverInterface::evalModel ---------------------------}" << std::endl;
#endif // HAVE_LIFEV_DEBUG
}

155
156
#endif /* HAVE_NOX_THYRA && HAVE_TRILINOS_RYTHMOS */

157
} // LifeV namespace