Albany: a Trilinos-based PDE code: Albany_ModelEvaluator.cpp Source File

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include "Albany_ModelEvaluator.hpp"
00008 #include "Teuchos_ScalarTraits.hpp"
00009 #include "Teuchos_TestForException.hpp"
00010 #include "Stokhos_EpetraVectorOrthogPoly.hpp"
00011 #include "Stokhos_EpetraMultiVectorOrthogPoly.hpp"
00012 #include "Stokhos_EpetraOperatorOrthogPoly.hpp"
00013 
00014 Albany::ModelEvaluator::ModelEvaluator(
00015   const Teuchos::RCP<Albany::Application>& app_,
00016   const Teuchos::RCP<Teuchos::ParameterList>& appParams) 
00017   : app(app_),
00018     supplies_prec(app_->suppliesPreconditioner())
00019 {
00020   Teuchos::RCP<Teuchos::FancyOStream> out = 
00021     Teuchos::VerboseObjectBase::getDefaultOStream();
00022 
00023   // Parameters (e.g., for sensitivities, SG expansions, ...)
00024   Teuchos::ParameterList& problemParams = appParams->sublist("Problem");
00025   Teuchos::ParameterList& parameterParams = 
00026     problemParams.sublist("Parameters");
00027   int num_param_vecs = 
00028     parameterParams.get("Number of Parameter Vectors", 0);
00029   bool using_old_parameter_list = false;
00030   if (parameterParams.isType<int>("Number")) {
00031     int numParameters = parameterParams.get<int>("Number");
00032     if (numParameters > 0) {
00033       num_param_vecs = 1;
00034       using_old_parameter_list = true;
00035     }
00036   }
00037   param_names.resize(num_param_vecs);
00038   *out << "Number of parameters vectors  = " << num_param_vecs << std::endl;
00039   for (int i=0; i<num_param_vecs; i++) {
00040     Teuchos::ParameterList* pList;
00041     if (using_old_parameter_list)
00042       pList = &parameterParams;
00043     else
00044       pList = &(parameterParams.sublist(Albany::strint("Parameter Vector",i)));
00045     int numParameters = pList->get<int>("Number");
00046     TEUCHOS_TEST_FOR_EXCEPTION(
00047       numParameters == 0, 
00048       Teuchos::Exceptions::InvalidParameter,
00049       std::endl << "Error!  In Albany::ModelEvaluator constructor:  " <<
00050       "Parameter vector " << i << " has zero parameters!" << std::endl);
00051     param_names[i] = 
00052       Teuchos::rcp(new Teuchos::Array<std::string>(numParameters));
00053     for (int j=0; j<numParameters; j++) {
00054       (*param_names[i])[j] = 
00055   pList->get<std::string>(Albany::strint("Parameter",j));
00056     }
00057     *out << "Number of parameters in parameter vector " << i << " = " 
00058    << numParameters << std::endl;
00059   }
00060 
00061   // Setup sacado and epetra storage for parameters  
00062   sacado_param_vec.resize(num_param_vecs);
00063   epetra_param_map.resize(num_param_vecs);
00064   epetra_param_vec.resize(num_param_vecs);
00065   p_sg_vals.resize(num_param_vecs);
00066   p_mp_vals.resize(num_param_vecs);
00067   const Epetra_Comm& comm = app->getMap()->Comm();
00068   for (int i=0; i<num_param_vecs; i++) {
00069 
00070     // Initialize Sacado parameter vector
00071     app->getParamLib()->fillVector<PHAL::AlbanyTraits::Residual>(
00072       *(param_names[i]), sacado_param_vec[i]);
00073 
00074     // Create Epetra map for parameter vector
00075     epetra_param_map[i] = 
00076       Teuchos::rcp(new Epetra_LocalMap((int) sacado_param_vec[i].size(), 0, comm));
00077 
00078     // Create Epetra vector for parameters
00079     epetra_param_vec[i] = 
00080       Teuchos::rcp(new Epetra_Vector(*(epetra_param_map[i])));
00081     for (unsigned int j=0; j<sacado_param_vec[i].size(); j++)
00082       (*(epetra_param_vec[i]))[j] = sacado_param_vec[i][j].baseValue;
00083 
00084     p_sg_vals[i].resize(sacado_param_vec[i].size());
00085     p_mp_vals[i].resize(sacado_param_vec[i].size());
00086   }  
00087 
00088   timer = Teuchos::TimeMonitor::getNewTimer("Albany: **Total Fill Time**");
00089 }
00090 
00091 Albany::ModelEvaluator::~ModelEvaluator(){
00092 #ifdef ALBANY_DEBUG
00093   std::cout << "Calling destructor for Albany_ModelEvaluator" << std::endl;
00094 #endif
00095 }
00096 
00097 // Overridden from EpetraExt::ModelEvaluator
00098 
00099 Teuchos::RCP<const Epetra_Map>
00100 Albany::ModelEvaluator::get_x_map() const
00101 {
00102   return app->getMap();
00103 }
00104 
00105 Teuchos::RCP<const Epetra_Map>
00106 Albany::ModelEvaluator::get_f_map() const
00107 {
00108   return app->getMap();
00109 }
00110 
00111 Teuchos::RCP<const Epetra_Map>
00112 Albany::ModelEvaluator::get_p_map(int l) const
00113 {
00114   TEUCHOS_TEST_FOR_EXCEPTION(
00115     l >= static_cast<int>(epetra_param_map.size()) || l < 0, 
00116     Teuchos::Exceptions::InvalidParameter,
00117     std::endl << 
00118     "Error!  Albany::ModelEvaluator::get_p_map():  " <<
00119     "Invalid parameter index l = " << l << std::endl);
00120 
00121   return epetra_param_map[l];
00122 }
00123 
00124 Teuchos::RCP<const Epetra_Map>
00125 Albany::ModelEvaluator::get_g_map(int l) const
00126 {
00127   TEUCHOS_TEST_FOR_EXCEPTION(
00128     l >= app->getNumResponses() || l < 0, 
00129     Teuchos::Exceptions::InvalidParameter,
00130     std::endl << 
00131     "Error!  Albany::ModelEvaluator::get_g_map():  " << 
00132     "Invalid response index l = " << l << std::endl);
00133 
00134   return app->getResponse(l)->responseMap();
00135 }
00136 
00137 Teuchos::RCP<const Teuchos::Array<std::string> >
00138 Albany::ModelEvaluator::get_p_names(int l) const
00139 {
00140   TEUCHOS_TEST_FOR_EXCEPTION(l >= static_cast<int>(param_names.size()) || l < 0, 
00141          Teuchos::Exceptions::InvalidParameter,
00142                      std::endl << 
00143                      "Error!  Albany::ModelEvaluator::get_p_names():  " <<
00144                      "Invalid parameter index l = " << l << std::endl);
00145 
00146   return param_names[l];
00147 }
00148 
00149 Teuchos::RCP<const Epetra_Vector>
00150 Albany::ModelEvaluator::get_x_init() const
00151 {
00152   return app->getAdaptSolMgr()->getInitialSolution();
00153 }
00154 
00155 Teuchos::RCP<const Epetra_Vector>
00156 Albany::ModelEvaluator::get_x_dot_init() const
00157 {
00158   return app->getAdaptSolMgr()->getInitialSolutionDot();
00159 }
00160 
00161 Teuchos::RCP<const Epetra_Vector>
00162 Albany::ModelEvaluator::get_x_dotdot_init() const
00163 {
00164   Teuchos::RCP<const Epetra_Vector> x_dotdot_init
00165     = Teuchos::rcp(new Epetra_Vector(app->getAdaptSolMgr()->getInitialSolutionDot()->Map()));
00166   
00167   return x_dotdot_init;
00168 }
00169 
00170 Teuchos::RCP<const Epetra_Vector>
00171 Albany::ModelEvaluator::get_p_init(int l) const
00172 {
00173   TEUCHOS_TEST_FOR_EXCEPTION(l >= static_cast<int>(param_names.size()) || l < 0, 
00174          Teuchos::Exceptions::InvalidParameter,
00175                      std::endl << 
00176                      "Error!  Albany::ModelEvaluator::get_p_init():  " <<
00177                      "Invalid parameter index l = " << l << std::endl);
00178   
00179   return epetra_param_vec[l];
00180 }
00181 
00182 Teuchos::RCP<Epetra_Operator>
00183 Albany::ModelEvaluator::create_W() const
00184 {
00185   return 
00186     Teuchos::rcp(new Epetra_CrsMatrix(::Copy, *(app->getJacobianGraph())));
00187 }
00188 
00189 Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner>
00190 Albany::ModelEvaluator::create_WPrec() const
00191 {
00192   Teuchos::RCP<Epetra_Operator> precOp = app->getPreconditioner();
00193 
00194   // Teko prec needs space for Jacobian as well
00195   Extra_W_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(create_W(), true);
00196 
00197   // bool is answer to: "Prec is already inverted?"
00198   return Teuchos::rcp(new EpetraExt::ModelEvaluator::Preconditioner(precOp,true));
00199 }
00200 
00201 Teuchos::RCP<Epetra_Operator>
00202 Albany::ModelEvaluator::create_DgDx_op(int j) const
00203 {
00204   TEUCHOS_TEST_FOR_EXCEPTION(
00205     j >= app->getNumResponses() || j < 0, 
00206     Teuchos::Exceptions::InvalidParameter,
00207     std::endl << 
00208     "Error!  Albany::ModelEvaluator::create_DgDx_op():  " << 
00209     "Invalid response index j = " << j << std::endl);
00210   
00211   return app->getResponse(j)->createGradientOp();
00212 }
00213 
00214 Teuchos::RCP<Epetra_Operator>
00215 Albany::ModelEvaluator::create_DgDx_dot_op(int j) const
00216 {
00217   TEUCHOS_TEST_FOR_EXCEPTION(
00218     j >= app->getNumResponses() || j < 0, 
00219     Teuchos::Exceptions::InvalidParameter,
00220     std::endl << 
00221     "Error!  Albany::ModelEvaluator::create_DgDx_dot_op():  " << 
00222     "Invalid response index j = " << j << std::endl);
00223   
00224   return app->getResponse(j)->createGradientOp();
00225 }
00226 
00227 Teuchos::RCP<Epetra_Operator>
00228 Albany::ModelEvaluator::create_DgDx_dotdot_op(int j) const
00229 {
00230   TEUCHOS_TEST_FOR_EXCEPTION(
00231     j >= app->getNumResponses() || j < 0, 
00232     Teuchos::Exceptions::InvalidParameter,
00233     std::endl << 
00234     "Error!  Albany::ModelEvaluator::create_DgDx_dotdot_op():  " << 
00235     "Invalid response index j = " << j << std::endl);
00236   
00237   return app->getResponse(j)->createGradientOp();
00238 }
00239 
00240 EpetraExt::ModelEvaluator::InArgs
00241 Albany::ModelEvaluator::createInArgs() const
00242 {
00243   InArgsSetup inArgs;
00244   inArgs.setModelEvalDescription(this->description());
00245 
00246   inArgs.setSupports(IN_ARG_t,true);
00247   inArgs.setSupports(IN_ARG_x,true);
00248   inArgs.setSupports(IN_ARG_x_dot,true);
00249   inArgs.setSupports(IN_ARG_x_dotdot,true);
00250   inArgs.setSupports(IN_ARG_alpha,true);
00251   inArgs.setSupports(IN_ARG_omega,true);
00252   inArgs.setSupports(IN_ARG_beta,true);
00253   inArgs.set_Np(param_names.size());
00254 
00255 #ifdef ALBANY_SG_MP
00256   inArgs.setSupports(IN_ARG_x_sg,true);
00257   inArgs.setSupports(IN_ARG_x_dot_sg,true);
00258   inArgs.setSupports(IN_ARG_x_dotdot_sg,true);
00259   for (int i=0; i<param_names.size(); i++)
00260     inArgs.setSupports(IN_ARG_p_sg, i, true);
00261   inArgs.setSupports(IN_ARG_sg_basis,true);
00262   inArgs.setSupports(IN_ARG_sg_quadrature,true);
00263   inArgs.setSupports(IN_ARG_sg_expansion,true);
00264   
00265   inArgs.setSupports(IN_ARG_x_mp,true);
00266   inArgs.setSupports(IN_ARG_x_dot_mp,true);
00267   inArgs.setSupports(IN_ARG_x_dotdot_mp,true);
00268   for (int i=0; i<param_names.size(); i++)
00269     inArgs.setSupports(IN_ARG_p_mp, i, true); 
00270 #endif
00271 
00272   return inArgs;
00273 }
00274 
00275 EpetraExt::ModelEvaluator::OutArgs
00276 Albany::ModelEvaluator::createOutArgs() const
00277 {
00278   OutArgsSetup outArgs;
00279   outArgs.setModelEvalDescription(this->description());
00280 
00281   int n_g = app->getNumResponses();
00282 
00283   // Deterministic
00284   outArgs.setSupports(OUT_ARG_f,true);
00285   outArgs.setSupports(OUT_ARG_W,true);
00286   outArgs.set_W_properties(
00287     DerivativeProperties(DERIV_LINEARITY_UNKNOWN, DERIV_RANK_FULL, true));
00288   if (supplies_prec) outArgs.setSupports(OUT_ARG_WPrec, true);
00289   outArgs.set_Np_Ng(param_names.size(), n_g);
00290 
00291   for (int i=0; i<param_names.size(); i++)
00292     outArgs.setSupports(OUT_ARG_DfDp, i, DerivativeSupport(DERIV_MV_BY_COL));
00293   for (int i=0; i<n_g; i++) {
00294     if (app->getResponse(i)->isScalarResponse()) {
00295       outArgs.setSupports(OUT_ARG_DgDx, i,
00296                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00297       outArgs.setSupports(OUT_ARG_DgDx_dot, i,
00298                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00299       outArgs.setSupports(OUT_ARG_DgDx_dotdot, i,
00300                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00301     }
00302     else {
00303       outArgs.setSupports(OUT_ARG_DgDx, i,
00304                           DerivativeSupport(DERIV_LINEAR_OP));
00305       outArgs.setSupports(OUT_ARG_DgDx_dot, i,
00306                           DerivativeSupport(DERIV_LINEAR_OP));
00307       outArgs.setSupports(OUT_ARG_DgDx_dotdot, i,
00308                           DerivativeSupport(DERIV_LINEAR_OP));
00309     }
00310 
00311     for (int j=0; j<param_names.size(); j++)
00312       outArgs.setSupports(OUT_ARG_DgDp, i, j,
00313                           DerivativeSupport(DERIV_MV_BY_COL));
00314   }
00315 
00316 
00317 #ifdef ALBANY_SG_MP
00318   // Stochastic
00319   outArgs.setSupports(OUT_ARG_f_sg,true);
00320   outArgs.setSupports(OUT_ARG_W_sg,true);
00321   for (int i=0; i<param_names.size(); i++)
00322     outArgs.setSupports(OUT_ARG_DfDp_sg, i, DerivativeSupport(DERIV_MV_BY_COL));
00323   for (int i=0; i<n_g; i++)
00324     outArgs.setSupports(OUT_ARG_g_sg, i, true);
00325   for (int i=0; i<n_g; i++) {
00326     if (app->getResponse(i)->isScalarResponse()) {
00327       outArgs.setSupports(OUT_ARG_DgDx_sg, i,
00328                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00329       outArgs.setSupports(OUT_ARG_DgDx_dot_sg, i,
00330                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00331       outArgs.setSupports(OUT_ARG_DgDx_dotdot_sg, i,
00332                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00333     }
00334     else {
00335       outArgs.setSupports(OUT_ARG_DgDx_sg, i,
00336                           DerivativeSupport(DERIV_LINEAR_OP));
00337       outArgs.setSupports(OUT_ARG_DgDx_dot_sg, i,
00338                           DerivativeSupport(DERIV_LINEAR_OP));
00339       outArgs.setSupports(OUT_ARG_DgDx_dotdot_sg, i,
00340                           DerivativeSupport(DERIV_LINEAR_OP));
00341     }
00342     for (int j=0; j<param_names.size(); j++)
00343       outArgs.setSupports(OUT_ARG_DgDp_sg, i, j,
00344                           DerivativeSupport(DERIV_MV_BY_COL));
00345   }
00346 
00347   // Multi-point
00348   outArgs.setSupports(OUT_ARG_f_mp,true);
00349   outArgs.setSupports(OUT_ARG_W_mp,true);
00350   for (int i=0; i<param_names.size(); i++)
00351     outArgs.setSupports(OUT_ARG_DfDp_mp, i, DerivativeSupport(DERIV_MV_BY_COL));
00352   for (int i=0; i<n_g; i++)
00353     outArgs.setSupports(OUT_ARG_g_mp, i, true);
00354   for (int i=0; i<n_g; i++) {
00355     outArgs.setSupports(OUT_ARG_g_mp, i, true);
00356     if (app->getResponse(i)->isScalarResponse()) {
00357       outArgs.setSupports(OUT_ARG_DgDx_mp, i,
00358                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00359       outArgs.setSupports(OUT_ARG_DgDx_dot_mp, i,
00360                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00361       outArgs.setSupports(OUT_ARG_DgDx_dotdot_mp, i,
00362                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00363     }
00364     else {
00365       outArgs.setSupports(OUT_ARG_DgDx_mp, i,
00366                           DerivativeSupport(DERIV_LINEAR_OP));
00367       outArgs.setSupports(OUT_ARG_DgDx_dot_mp, i,
00368                           DerivativeSupport(DERIV_LINEAR_OP));
00369       outArgs.setSupports(OUT_ARG_DgDx_dotdot_mp, i,
00370                           DerivativeSupport(DERIV_LINEAR_OP));
00371     }
00372     for (int j=0; j<param_names.size(); j++)
00373       outArgs.setSupports(OUT_ARG_DgDp_mp, i, j,
00374                           DerivativeSupport(DERIV_MV_BY_COL));
00375   }
00376 #endif
00377 
00378   return outArgs;
00379 }
00380 
00381 void 
00382 Albany::ModelEvaluator::evalModel(const InArgs& inArgs, 
00383          const OutArgs& outArgs) const
00384 {
00385   Teuchos::TimeMonitor Timer(*timer); //start timer
00386   //
00387   // Get the input arguments
00388   //
00389   Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
00390   Teuchos::RCP<const Epetra_Vector> x_dot;
00391   Teuchos::RCP<const Epetra_Vector> x_dotdot;
00392   double alpha     = 0.0;
00393   double omega     = 0.0;
00394   double beta      = 1.0;
00395   double curr_time = 0.0;
00396   x_dot = inArgs.get_x_dot();
00397   x_dotdot = inArgs.get_x_dotdot();
00398   if (x_dot != Teuchos::null || x_dotdot != Teuchos::null) {
00399     alpha = inArgs.get_alpha();
00400     omega = inArgs.get_omega();
00401     beta = inArgs.get_beta();
00402     curr_time  = inArgs.get_t();
00403   }
00404   for (int i=0; i<inArgs.Np(); i++) {
00405     Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i);
00406     if (p != Teuchos::null) {
00407       for (unsigned int j=0; j<sacado_param_vec[i].size(); j++)
00408   sacado_param_vec[i][j].baseValue = (*p)[j];
00409     }
00410   }
00411 
00412   //
00413   // Get the output arguments
00414   //
00415   EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out = outArgs.get_f();
00416   Teuchos::RCP<Epetra_Operator> W_out = outArgs.get_W();
00417 
00418   // Cast W to a CrsMatrix, throw an exception if this fails
00419   Teuchos::RCP<Epetra_CrsMatrix> W_out_crs;
00420 
00421   if (W_out != Teuchos::null)
00422     W_out_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out, true);
00423 
00424 int test_var = 0;
00425 if(test_var != 0){
00426 std::cout << "The current solution length is: " << x->MyLength() << std::endl;
00427 x->Print(std::cout);
00428 
00429 }
00430   
00431   // Get preconditioner operator, if requested
00432   Teuchos::RCP<Epetra_Operator> WPrec_out;
00433   if (outArgs.supports(OUT_ARG_WPrec)) WPrec_out = outArgs.get_WPrec();
00434 
00435   //
00436   // Compute the functions
00437   //
00438   bool f_already_computed = false;
00439 
00440   // W matrix
00441   if (W_out != Teuchos::null) {
00442     app->computeGlobalJacobian(alpha, beta, omega, curr_time, x_dot.get(), x_dotdot.get(),*x, 
00443              sacado_param_vec, f_out.get(), *W_out_crs);
00444     f_already_computed=true;
00445 if(test_var != 0){
00446 //std::cout << "The current rhs length is: " << f_out->MyLength() << std::endl;
00447 //f_out->Print(std::cout);
00448 std::cout << "The current Jacobian length is: " << W_out_crs->NumGlobalRows() << std::endl;
00449 W_out_crs->Print(std::cout);
00450 }
00451   }
00452 
00453   if (WPrec_out != Teuchos::null) {
00454     app->computeGlobalJacobian(alpha, beta, omega, curr_time, x_dot.get(), x_dotdot.get(), *x, 
00455              sacado_param_vec, f_out.get(), *Extra_W_crs);
00456     f_already_computed=true;
00457 if(test_var != 0){
00458 //std::cout << "The current rhs length is: " << f_out->MyLength() << std::endl;
00459 //f_out->Print(std::cout);
00460 std::cout << "The current preconditioner length is: " << Extra_W_crs->NumGlobalRows() << std::endl;
00461 Extra_W_crs->Print(std::cout);
00462 }
00463 
00464     app->computeGlobalPreconditioner(Extra_W_crs, WPrec_out);
00465   }
00466 
00467   // df/dp
00468   for (int i=0; i<outArgs.Np(); i++) {
00469     Teuchos::RCP<Epetra_MultiVector> dfdp_out = 
00470       outArgs.get_DfDp(i).getMultiVector();
00471     if (dfdp_out != Teuchos::null) {
00472       Teuchos::Array<int> p_indexes = 
00473         outArgs.get_DfDp(i).getDerivativeMultiVector().getParamIndexes();
00474       Teuchos::RCP<ParamVec> p_vec;
00475       if (p_indexes.size() == 0)
00476         p_vec = Teuchos::rcp(&sacado_param_vec[i],false);
00477       else {
00478         p_vec = Teuchos::rcp(new ParamVec);
00479         for (int j=0; j<p_indexes.size(); j++)
00480           p_vec->addParam(sacado_param_vec[i][p_indexes[j]].family, 
00481                           sacado_param_vec[i][p_indexes[j]].baseValue);
00482       }
00483 
00484       app->computeGlobalTangent(0.0, 0.0, 0.0, curr_time, false, x_dot.get(), x_dotdot.get(), *x, 
00485                                 sacado_param_vec, p_vec.get(),
00486                                 NULL, NULL, NULL, NULL, f_out.get(), NULL, 
00487                                 dfdp_out.get());
00488 
00489       f_already_computed=true;
00490 if(test_var != 0){
00491 std::cout << "The current rhs length is: " << f_out->MyLength() << std::endl;
00492 f_out->Print(std::cout);
00493 }
00494     }
00495   }
00496 
00497   // f
00498   if (app->is_adjoint) {
00499     Derivative f_deriv(f_out, DERIV_TRANS_MV_BY_ROW);
00500     int response_index = 0; // need to add capability for sending this in
00501     app->evaluateResponseDerivative(response_index, curr_time, x_dot.get(), x_dotdot.get(), *x, 
00502             sacado_param_vec, NULL, 
00503             NULL, f_deriv, Derivative(), Derivative(), Derivative());
00504   }
00505   else {
00506     if (f_out != Teuchos::null && !f_already_computed) {
00507       app->computeGlobalResidual(curr_time, x_dot.get(), x_dotdot.get(), *x, 
00508                  sacado_param_vec, *f_out);
00509 if(test_var != 0){
00510 std::cout << "The current rhs length is: " << f_out->MyLength() << std::endl;
00511 f_out->Print(std::cout);
00512 }
00513     }
00514   }
00515 
00516   // Response functions
00517   for (int i=0; i<outArgs.Ng(); i++) {
00518     Teuchos::RCP<Epetra_Vector> g_out = outArgs.get_g(i);
00519     bool g_computed = false;
00520 
00521     Derivative dgdx_out = outArgs.get_DgDx(i);
00522     Derivative dgdxdot_out = outArgs.get_DgDx_dot(i);
00523     Derivative dgdxdotdot_out = outArgs.get_DgDx_dotdot(i);
00524 
00525     // dg/dx, dg/dxdot
00526     if (!dgdx_out.isEmpty() || !dgdxdot_out.isEmpty() || !dgdxdotdot_out.isEmpty() ) {
00527       app->evaluateResponseDerivative(i, curr_time, x_dot.get(), x_dotdot.get(), *x,
00528                                       sacado_param_vec, NULL,
00529                                       g_out.get(), dgdx_out,
00530                                       dgdxdot_out, dgdxdotdot_out, Derivative());
00531       g_computed = true;
00532     }
00533 
00534     // dg/dp
00535     for (int j=0; j<outArgs.Np(); j++) {
00536       Teuchos::RCP<Epetra_MultiVector> dgdp_out =
00537         outArgs.get_DgDp(i,j).getMultiVector();
00538       if (dgdp_out != Teuchos::null) {
00539         Teuchos::Array<int> p_indexes =
00540           outArgs.get_DgDp(i,j).getDerivativeMultiVector().getParamIndexes();
00541         Teuchos::RCP<ParamVec> p_vec;
00542         if (p_indexes.size() == 0)
00543           p_vec = Teuchos::rcp(&sacado_param_vec[j],false);
00544         else {
00545           p_vec = Teuchos::rcp(new ParamVec);
00546           for (int k=0; k<p_indexes.size(); k++)
00547             p_vec->addParam(sacado_param_vec[j][p_indexes[k]].family,
00548                             sacado_param_vec[j][p_indexes[k]].baseValue);
00549         }
00550         app->evaluateResponseTangent(i, alpha, beta, omega, curr_time, false,
00551                                      x_dot.get(), x_dotdot.get(), *x,
00552                                      sacado_param_vec, p_vec.get(),
00553                                      NULL, NULL, NULL, NULL, g_out.get(), NULL,
00554                                      dgdp_out.get());
00555         g_computed = true;
00556       }
00557     }
00558 
00559     if (g_out != Teuchos::null && !g_computed)
00560       app->evaluateResponse(i, curr_time, x_dot.get(), x_dotdot.get(), *x, sacado_param_vec, 
00561           *g_out);
00562   }
00563 
00564   //
00565   // Stochastic Galerkin
00566   //
00567 #ifdef ALBANY_SG_MP
00568   InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
00569   if (x_sg != Teuchos::null) {
00570     app->init_sg(inArgs.get_sg_basis(), 
00571      inArgs.get_sg_quadrature(), 
00572      inArgs.get_sg_expansion(), 
00573      x_sg->productComm());
00574     InArgs::sg_const_vector_t x_dot_sg  = inArgs.get_x_dot_sg();
00575     InArgs::sg_const_vector_t x_dotdot_sg  = inArgs.get_x_dotdot_sg();
00576     InArgs::sg_const_vector_t epetra_p_sg = inArgs.get_p_sg(0);
00577     Teuchos::Array<int> p_sg_index;
00578     for (int i=0; i<inArgs.Np(); i++) {
00579       InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(i);
00580       if (p_sg != Teuchos::null) {
00581   p_sg_index.push_back(i);
00582   for (int j=0; j<p_sg_vals[i].size(); j++) {
00583     int num_sg_blocks = p_sg->size();
00584     p_sg_vals[i][j].reset(app->getStochasticExpansion(), num_sg_blocks);
00585     p_sg_vals[i][j].copyForWrite();
00586     for (int l=0; l<num_sg_blocks; l++) {
00587       p_sg_vals[i][j].fastAccessCoeff(l) = (*p_sg)[l][j];
00588     }
00589   }
00590       }
00591     }
00592 
00593     OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
00594     OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
00595     bool f_sg_computed = false;
00596 
00597     // W_sg
00598     if (W_sg != Teuchos::null) {
00599       Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> W_sg_crs(W_sg->basis(), 
00600                  W_sg->map());
00601       for (int i=0; i<W_sg->size(); i++)
00602   W_sg_crs.setCoeffPtr(
00603     i,
00604     Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_sg->getCoeffPtr(i)));
00605       app->computeGlobalSGJacobian(alpha, beta, omega, curr_time, 
00606            x_dot_sg.get(),  x_dotdot_sg.get(), *x_sg, 
00607            sacado_param_vec, p_sg_index, p_sg_vals,
00608            f_sg.get(), W_sg_crs);
00609       f_sg_computed = true;
00610     }
00611 
00612     // df/dp_sg
00613     for (int i=0; i<outArgs.Np(); i++) {
00614       Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > dfdp_sg
00615         = outArgs.get_DfDp_sg(i).getMultiVector();
00616       if (dfdp_sg != Teuchos::null) {
00617         Teuchos::Array<int> p_indexes =
00618           outArgs.get_DfDp_sg(i).getDerivativeMultiVector().getParamIndexes();
00619         Teuchos::RCP<ParamVec> p_vec;
00620         if (p_indexes.size() == 0)
00621           p_vec = Teuchos::rcp(&sacado_param_vec[i],false);
00622         else {
00623           p_vec = Teuchos::rcp(new ParamVec);
00624           for (int j=0; j<p_indexes.size(); j++)
00625             p_vec->addParam(sacado_param_vec[i][p_indexes[j]].family,
00626                             sacado_param_vec[i][p_indexes[j]].baseValue);
00627         }
00628 
00629         app->computeGlobalSGTangent(0.0, 0.0, 0.0, curr_time, false,
00630                                     x_dot_sg.get(), x_dotdot_sg.get(),*x_sg,
00631                                     sacado_param_vec, p_sg_index, p_sg_vals,
00632                                     p_vec.get(), NULL, NULL, NULL, NULL,
00633                                     f_sg.get(), NULL, dfdp_sg.get());
00634 
00635         f_sg_computed = true;
00636       }
00637     }
00638 
00639     if (f_sg != Teuchos::null && !f_sg_computed)
00640       app->computeGlobalSGResidual(curr_time, x_dot_sg.get(), x_dotdot_sg.get(),*x_sg, 
00641            sacado_param_vec, p_sg_index, p_sg_vals,
00642            *f_sg);
00643 
00644     // Response functions
00645     for (int i=0; i<outArgs.Ng(); i++) {
00646       OutArgs::sg_vector_t g_sg = outArgs.get_g_sg(i);
00647       bool g_sg_computed = false;
00648 
00649       SGDerivative dgdx_sg = outArgs.get_DgDx_sg(i);
00650       SGDerivative dgdxdot_sg = outArgs.get_DgDx_dot_sg(i);
00651       SGDerivative dgdxdotdot_sg = outArgs.get_DgDx_dotdot_sg(i);
00652 
00653       // dg/dx, dg/dxdot
00654       if (!dgdx_sg.isEmpty() || !dgdxdot_sg.isEmpty() || !dgdxdotdot_sg.isEmpty()) {
00655         app->evaluateSGResponseDerivative(
00656             i, curr_time, x_dot_sg.get(), x_dotdot_sg.get(), *x_sg,
00657             sacado_param_vec, p_sg_index, p_sg_vals,
00658             NULL, g_sg.get(), dgdx_sg,
00659             dgdxdot_sg, dgdxdotdot_sg, SGDerivative());
00660         g_sg_computed = true;
00661       }
00662 
00663       // dg/dp
00664       for (int j=0; j<outArgs.Np(); j++) {
00665         Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > dgdp_sg =
00666           outArgs.get_DgDp_sg(i,j).getMultiVector();
00667         if (dgdp_sg != Teuchos::null) {
00668           Teuchos::Array<int> p_indexes =
00669             outArgs.get_DgDp_sg(i,j).getDerivativeMultiVector().getParamIndexes();
00670           Teuchos::RCP<ParamVec> p_vec;
00671           if (p_indexes.size() == 0)
00672             p_vec = Teuchos::rcp(&sacado_param_vec[j],false);
00673           else {
00674             p_vec = Teuchos::rcp(new ParamVec);
00675             for (int k=0; k<p_indexes.size(); k++)
00676               p_vec->addParam(sacado_param_vec[j][p_indexes[k]].family,
00677                               sacado_param_vec[j][p_indexes[k]].baseValue);
00678           }
00679           app->evaluateSGResponseTangent(i, alpha, beta, omega, curr_time, false,
00680                                          x_dot_sg.get(), x_dotdot_sg.get(), *x_sg,
00681                                          sacado_param_vec, p_sg_index,
00682                                          p_sg_vals, p_vec.get(),
00683                                          NULL, NULL, NULL, NULL, g_sg.get(),
00684                                          NULL, dgdp_sg.get());
00685           g_sg_computed = true;
00686 
00687         }
00688       }
00689 
00690       if (g_sg != Teuchos::null && !g_sg_computed)
00691   app->evaluateSGResponse(i, curr_time, x_dot_sg.get(), x_dotdot_sg.get(), *x_sg, 
00692         sacado_param_vec, p_sg_index, p_sg_vals, 
00693         *g_sg);
00694     }
00695   }
00696 
00697   //
00698   // Multi-point evaluation
00699   //
00700   mp_const_vector_t x_mp = inArgs.get_x_mp();
00701   if (x_mp != Teuchos::null) {
00702     mp_const_vector_t x_dot_mp  = inArgs.get_x_dot_mp();
00703     mp_const_vector_t x_dotdot_mp  = inArgs.get_x_dotdot_mp();
00704     Teuchos::Array<int> p_mp_index;
00705     for (int i=0; i<inArgs.Np(); i++) {
00706       mp_const_vector_t p_mp = inArgs.get_p_mp(i);
00707       if (p_mp != Teuchos::null) {
00708   p_mp_index.push_back(i);
00709   for (int j=0; j<p_mp_vals[i].size(); j++) {
00710     int num_mp_blocks = p_mp->size();
00711     p_mp_vals[i][j].reset(num_mp_blocks);
00712     p_mp_vals[i][j].copyForWrite();
00713     for (int l=0; l<num_mp_blocks; l++) {
00714       p_mp_vals[i][j].fastAccessCoeff(l) = (*p_mp)[l][j];
00715     }
00716   }
00717       }
00718     }
00719     
00720     mp_vector_t f_mp = outArgs.get_f_mp();
00721     mp_operator_t W_mp = outArgs.get_W_mp();
00722     bool f_mp_computed = false;
00723 
00724     // W_mp
00725     if (W_mp != Teuchos::null) {
00726       Stokhos::ProductContainer<Epetra_CrsMatrix> W_mp_crs(W_mp->map());
00727       for (int i=0; i<W_mp->size(); i++)
00728   W_mp_crs.setCoeffPtr(
00729     i,
00730     Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_mp->getCoeffPtr(i)));
00731       app->computeGlobalMPJacobian(alpha, beta, omega, curr_time, 
00732            x_dot_mp.get(), x_dotdot_mp.get(), *x_mp, 
00733            sacado_param_vec, p_mp_index, p_mp_vals,
00734            f_mp.get(), W_mp_crs);
00735       f_mp_computed = true;
00736     }
00737 
00738     // df/dp_mp
00739     for (int i=0; i<outArgs.Np(); i++) {
00740       Teuchos::RCP< Stokhos::ProductEpetraMultiVector > dfdp_mp
00741         = outArgs.get_DfDp_mp(i).getMultiVector();
00742       if (dfdp_mp != Teuchos::null) {
00743         Teuchos::Array<int> p_indexes =
00744           outArgs.get_DfDp_mp(i).getDerivativeMultiVector().getParamIndexes();
00745         Teuchos::RCP<ParamVec> p_vec;
00746         if (p_indexes.size() == 0)
00747           p_vec = Teuchos::rcp(&sacado_param_vec[i],false);
00748         else {
00749           p_vec = Teuchos::rcp(new ParamVec);
00750           for (int j=0; j<p_indexes.size(); j++)
00751             p_vec->addParam(sacado_param_vec[i][p_indexes[j]].family,
00752                             sacado_param_vec[i][p_indexes[j]].baseValue);
00753         }
00754 
00755         app->computeGlobalMPTangent(0.0, 0.0, 0.0, curr_time, false,
00756                                     x_dot_mp.get(), x_dotdot_mp.get(), *x_mp,
00757                                     sacado_param_vec, p_mp_index, p_mp_vals,
00758                                     p_vec.get(), NULL, NULL, NULL, NULL,
00759                                     f_mp.get(), NULL, dfdp_mp.get());
00760 
00761         f_mp_computed = true;
00762       }
00763     }
00764 
00765     if (f_mp != Teuchos::null && !f_mp_computed)
00766       app->computeGlobalMPResidual(curr_time, x_dot_mp.get(), x_dotdot_mp.get(), *x_mp, 
00767            sacado_param_vec, p_mp_index, p_mp_vals,
00768            *f_mp);
00769 
00770     // Response functions
00771     for (int i=0; i<outArgs.Ng(); i++) {
00772       mp_vector_t g_mp = outArgs.get_g_mp(i);
00773       bool g_mp_computed = false;
00774 
00775       MPDerivative dgdx_mp = outArgs.get_DgDx_mp(i);
00776       MPDerivative dgdxdot_mp = outArgs.get_DgDx_dot_mp(i);
00777       MPDerivative dgdxdotdot_mp = outArgs.get_DgDx_dotdot_mp(i);
00778 
00779       // dg/dx, dg/dxdot
00780       if (!dgdx_mp.isEmpty() || !dgdxdot_mp.isEmpty() || !dgdxdotdot_mp.isEmpty() ) {
00781         app->evaluateMPResponseDerivative(
00782             i, curr_time, x_dot_mp.get(), x_dotdot_mp.get(), *x_mp,
00783             sacado_param_vec, p_mp_index, p_mp_vals,
00784             NULL, g_mp.get(), dgdx_mp,
00785             dgdxdot_mp, dgdxdotdot_mp, MPDerivative());
00786         g_mp_computed = true;
00787       }
00788 
00789       // dg/dp
00790       for (int j=0; j<outArgs.Np(); j++) {
00791         Teuchos::RCP< Stokhos::ProductEpetraMultiVector > dgdp_mp =
00792           outArgs.get_DgDp_mp(i,j).getMultiVector();
00793         if (dgdp_mp != Teuchos::null) {
00794           Teuchos::Array<int> p_indexes =
00795             outArgs.get_DgDp_mp(i,j).getDerivativeMultiVector().getParamIndexes();
00796           Teuchos::RCP<ParamVec> p_vec;
00797           if (p_indexes.size() == 0)
00798             p_vec = Teuchos::rcp(&sacado_param_vec[j],false);
00799           else {
00800             p_vec = Teuchos::rcp(new ParamVec);
00801             for (int k=0; k<p_indexes.size(); k++)
00802               p_vec->addParam(sacado_param_vec[j][p_indexes[k]].family,
00803                               sacado_param_vec[j][p_indexes[k]].baseValue);
00804           }
00805           app->evaluateMPResponseTangent(i, alpha, beta, omega, curr_time, false,
00806                                          x_dot_mp.get(), x_dotdot_mp.get(), *x_mp,
00807                                          sacado_param_vec, p_mp_index,
00808                                          p_mp_vals, p_vec.get(),
00809                                          NULL, NULL, NULL, NULL, g_mp.get(),
00810                                          NULL, dgdp_mp.get());
00811           g_mp_computed = true;
00812         }
00813       }
00814 
00815       if (g_mp != Teuchos::null && !g_mp_computed)
00816   app->evaluateMPResponse(i, curr_time, x_dot_mp.get(), x_dotdot_mp.get(), *x_mp, 
00817         sacado_param_vec, p_mp_index, p_mp_vals, 
00818         *g_mp);
00819     }
00820   }
00821 #endif //ALBANY_SG_MP
00822 }