00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "LOCA_Epetra_ModelEvaluatorInterface.H"
00043 #include "LOCA_Parameter_Vector.H"
00044 #include "NOX_Epetra_MultiVector.H"
00045 #include "LOCA_MultiContinuation_AbstractGroup.H"
00046
00047 #include "EpetraExt_ModelEvaluator.h"
00048 #include "Epetra_Operator.h"
00049
00050
00051
00052 LOCA::Epetra::ModelEvaluatorInterface::
00053 ModelEvaluatorInterface(
00054 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00055 const Teuchos::RefCountPtr<EpetraExt::ModelEvaluator>& m,
00056 double perturb) :
00057 NOX::Epetra::ModelEvaluatorInterface(m),
00058 LOCA::DerivUtils(global_data, perturb),
00059 param_vec(*(m->get_p_init(0))),
00060 loca_param_vec(),
00061 x_dot(NULL)
00062 {
00063
00064 Teuchos::RefCountPtr<const Teuchos::Array<std::string> > param_names =
00065 m->get_p_names(0);
00066 for (std::size_t i=0; i<param_names->size(); i++)
00067 loca_param_vec.addParameter((*param_names)[i], param_vec[i]);
00068 }
00069
00070
00071
00072 LOCA::Epetra::ModelEvaluatorInterface::~ModelEvaluatorInterface()
00073 {
00074 if (x_dot)
00075 delete x_dot;
00076 }
00077
00078
00079
00080 const LOCA::ParameterVector&
00081 LOCA::Epetra::ModelEvaluatorInterface::getLOCAParameterVector() const
00082 {
00083 return loca_param_vec;
00084 }
00085
00086
00087
00088 bool LOCA::Epetra::ModelEvaluatorInterface::
00089 computeF(const Epetra_Vector& x, Epetra_Vector& F, const FillType fillFlag)
00090 {
00091
00092 EpetraExt::ModelEvaluator::InArgs inargs = model_->createInArgs();
00093 inargs.set_x(Teuchos::rcp(&x, false));
00094 inargs.set_p(0, Teuchos::rcp(¶m_vec, false));
00095 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_x_dot)) {
00096
00097 if (x_dot == NULL)
00098 x_dot = new Epetra_Vector(x.Map());
00099 inargs.set_x_dot(Teuchos::rcp(x_dot, false));
00100 }
00101
00102
00103 EpetraExt::ModelEvaluator::OutArgs outargs = model_->createOutArgs();
00104 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> eval_f;
00105 Teuchos::RefCountPtr<Epetra_Vector> f = Teuchos::rcp(&F, false);
00106 if (fillFlag == NOX::Epetra::Interface::Required::Residual)
00107 eval_f.reset(f, EpetraExt::ModelEvaluator::EVAL_TYPE_EXACT);
00108 else if (fillFlag == NOX::Epetra::Interface::Required::Jac)
00109 eval_f.reset(f, EpetraExt::ModelEvaluator::EVAL_TYPE_APPROX_DERIV);
00110 else
00111 eval_f.reset(f, EpetraExt::ModelEvaluator::EVAL_TYPE_VERY_APPROX_DERIV);
00112 outargs.set_f(eval_f);
00113
00114 model_->evalModel(inargs, outargs);
00115
00116 return true;
00117 }
00118
00119
00120
00121 bool LOCA::Epetra::ModelEvaluatorInterface::
00122 computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)
00123 {
00124
00125 EpetraExt::ModelEvaluator::InArgs inargs = model_->createInArgs();
00126 inargs.set_x(Teuchos::rcp(&x, false));
00127 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_alpha))
00128 inargs.set_alpha(0.0);
00129 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_beta))
00130 inargs.set_beta(1.0);
00131 inargs.set_p(0, Teuchos::rcp(¶m_vec, false));
00132 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_x_dot)) {
00133
00134 if (x_dot == NULL)
00135 x_dot = new Epetra_Vector(x.Map());
00136 inargs.set_x_dot(Teuchos::rcp(x_dot, false));
00137 }
00138
00139
00140 EpetraExt::ModelEvaluator::OutArgs outargs = model_->createOutArgs();
00141 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> eval_f;
00142 outargs.set_f(eval_f);
00143 outargs.set_W(Teuchos::rcp(&Jac, false));
00144
00145 model_->evalModel(inargs, outargs);
00146
00147 return true;
00148 }
00149
00150
00151
00152 bool LOCA::Epetra::ModelEvaluatorInterface::
00153 computePreconditioner(const Epetra_Vector& x,
00154 Epetra_Operator& M,
00155 Teuchos::ParameterList* precParams)
00156 {
00157
00158 EpetraExt::ModelEvaluator::InArgs inargs = model_->createInArgs();
00159 inargs.set_x(Teuchos::rcp(&x, false));
00160 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_alpha))
00161 inargs.set_alpha(0.0);
00162 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_beta))
00163 inargs.set_beta(1.0);
00164 inargs.set_p(0, Teuchos::rcp(¶m_vec, false));
00165 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_x_dot)) {
00166
00167 if (x_dot == NULL)
00168 x_dot = new Epetra_Vector(x.Map());
00169 inargs.set_x_dot(Teuchos::rcp(x_dot, false));
00170 }
00171
00172
00173 EpetraExt::ModelEvaluator::OutArgs outargs = model_->createOutArgs();
00174 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> eval_f;
00175 outargs.set_f(eval_f);
00176 outargs.set_W(Teuchos::rcp(&M, false));
00177
00178 model_->evalModel(inargs, outargs);
00179
00180 return true;
00181 }
00182
00183
00184
00185 void LOCA::Epetra::ModelEvaluatorInterface::
00186 setParameters(const ParameterVector& p)
00187 {
00188 for (int i=0; i<p.length(); i++)
00189 param_vec[i] = p[i];
00190 }
00191
00192
00193
00194 bool LOCA::Epetra::ModelEvaluatorInterface::
00195 computeShiftedMatrix(double alpha, double beta, const Epetra_Vector& x,
00196 Epetra_Operator& A)
00197 {
00198
00199 EpetraExt::ModelEvaluator::InArgs inargs = model_->createInArgs();
00200 inargs.set_x(Teuchos::rcp(&x, false));
00201 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_alpha))
00202 inargs.set_alpha(-beta);
00203 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_beta))
00204 inargs.set_beta(alpha);
00205 inargs.set_p(0, Teuchos::rcp(¶m_vec, false));
00206 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_x_dot)) {
00207
00208 if (x_dot == NULL)
00209 x_dot = new Epetra_Vector(x.Map());
00210 inargs.set_x_dot(Teuchos::rcp(x_dot, false));
00211 }
00212
00213
00214 EpetraExt::ModelEvaluator::OutArgs outargs = model_->createOutArgs();
00215 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> eval_f;
00216 outargs.set_f(eval_f);
00217 outargs.set_W(Teuchos::rcp(&A, false));
00218
00219 model_->evalModel(inargs, outargs);
00220
00221 return true;
00222 }
00223
00224 LOCA::Epetra::ModelEvaluatorInterface::
00225 ModelEvaluatorInterface(const LOCA::Epetra::ModelEvaluatorInterface& m) :
00226 NOX::Epetra::ModelEvaluatorInterface(m),
00227 LOCA::DerivUtils(m),
00228 param_vec(m.param_vec),
00229 loca_param_vec(m.loca_param_vec),
00230 x_dot(NULL)
00231 {
00232 if (m.x_dot != NULL) {
00233 x_dot = new Epetra_Vector(*m.x_dot);
00234 }
00235 }
00236
00237 Teuchos::RCP<LOCA::DerivUtils>
00238 LOCA::Epetra::ModelEvaluatorInterface::
00239 clone(NOX::CopyType type) const
00240 {
00241 return Teuchos::rcp(new LOCA::Epetra::ModelEvaluatorInterface(*this));
00242 }
00243
00244 NOX::Abstract::Group::ReturnType
00245 LOCA::Epetra::ModelEvaluatorInterface::
00246 computeDfDp(LOCA::MultiContinuation::AbstractGroup& grp,
00247 const vector<int>& param_ids,
00248 NOX::Abstract::MultiVector& result,
00249 bool isValidF) const
00250 {
00251
00252 NOX::Epetra::Vector& f = dynamic_cast<NOX::Epetra::Vector&>(result[0]);
00253 Epetra_Vector& epetra_f = f.getEpetraVector();
00254
00255 std::vector<int> dfdp_index(result.numVectors()-1);
00256 for (unsigned int i=0; i<dfdp_index.size(); i++)
00257 dfdp_index[i] = i+1;
00258 Teuchos::RefCountPtr<NOX::Epetra::MultiVector> dfdp =
00259 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(result.subView(dfdp_index));
00260 Epetra_MultiVector& epetra_dfdp = dfdp->getEpetraMultiVector();
00261
00262
00263 EpetraExt::ModelEvaluator::InArgs inargs = model_->createInArgs();
00264 const NOX::Epetra::Vector& x =
00265 dynamic_cast<const NOX::Epetra::Vector&>(grp.getX());
00266 const Epetra_Vector& epetra_x = x.getEpetraVector();
00267 inargs.set_x(Teuchos::rcp(&epetra_x, false));
00268 inargs.set_p(0, Teuchos::rcp(¶m_vec, false));
00269 if (inargs.supports(EpetraExt::ModelEvaluator::IN_ARG_x_dot)) {
00270
00271 if (x_dot == NULL)
00272 x_dot = new Epetra_Vector(epetra_x.Map());
00273 inargs.set_x_dot(Teuchos::rcp(x_dot, false));
00274 }
00275
00276
00277 EpetraExt::ModelEvaluator::OutArgs outargs = model_->createOutArgs();
00278 if (!isValidF) {
00279 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> eval_f;
00280 Teuchos::RefCountPtr<Epetra_Vector> F = Teuchos::rcp(&epetra_f, false);
00281 eval_f.reset(F, EpetraExt::ModelEvaluator::EVAL_TYPE_EXACT);
00282 outargs.set_f(eval_f);
00283 }
00284 Teuchos::RefCountPtr<Epetra_MultiVector> DfDp =
00285 Teuchos::rcp(&epetra_dfdp, false);
00286 Teuchos::Array<int> param_indexes(param_ids.size());
00287 for (unsigned int i=0; i<param_ids.size(); i++)
00288 param_indexes[i] = param_ids[i];
00289 EpetraExt::ModelEvaluator::DerivativeMultiVector dmv(DfDp, EpetraExt::ModelEvaluator::DERIV_MV_BY_COL,
00290 param_indexes);
00291 EpetraExt::ModelEvaluator::Derivative deriv(dmv);
00292 outargs.set_DfDp(0, deriv);
00293
00294 model_->evalModel(inargs, outargs);
00295
00296 return NOX::Abstract::Group::Ok;
00297 }