FEApp_ModelEvaluator.cpp

Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 // @HEADER
00031 
00032 #include "FEApp_ModelEvaluator.hpp"
00033 #include "Teuchos_ScalarTraits.hpp"
00034 #include "Teuchos_TestForException.hpp"
00035 
00036 FEApp::ModelEvaluator::ModelEvaluator(
00037   const Teuchos::RCP<FEApp::Application>& app_,
00038   const Teuchos::RCP< Teuchos::Array<std::string> >& free_param_names,
00039   const Teuchos::RCP< Teuchos::Array<std::string> >& sg_param_names) 
00040   : app(app_),
00041     supports_p(false),
00042     supports_g(false),
00043     supports_sg(false),
00044     eval_W_with_f(false)
00045 {
00046   // Compute number of parameter vectors
00047   int num_param_vecs = 0;
00048   if (free_param_names != Teuchos::null) {
00049     supports_p = true;
00050     num_param_vecs = 1;
00051   }
00052 #if SG_ACTIVE
00053   if (sg_param_names != Teuchos::null) {
00054     supports_sg = true;
00055     num_param_vecs = 2;
00056   }
00057 #endif
00058 
00059   if (num_param_vecs > 0) {
00060     param_names.resize(num_param_vecs);
00061     sacado_param_vec.resize(num_param_vecs);
00062     epetra_param_map.resize(num_param_vecs);
00063     epetra_param_vec.resize(num_param_vecs);
00064 
00065     // Set parameter names
00066     param_names[0] = free_param_names;
00067     if (num_param_vecs == 2)
00068       param_names[1] = sg_param_names;
00069 
00070     // Initialize each parameter vector
00071     const Epetra_Comm& comm = app->getMap()->Comm();
00072     for (int i=0; i<num_param_vecs; i++) {
00073 
00074       // Initialize Sacado parameter vector
00075       sacado_param_vec[i] = Teuchos::rcp(new ParamVec);
00076       if (param_names[i] != Teuchos::null)
00077   app->getParamLib()->fillVector<FEApp::ResidualType>(*(param_names[i]), 
00078                   *(sacado_param_vec[i]));
00079 
00080       // Create Epetra map for parameter vector
00081       epetra_param_map[i] = 
00082   Teuchos::rcp(new Epetra_LocalMap(sacado_param_vec[i]->size(), 0, comm));
00083 
00084       // Create Epetra vector for parameters
00085       epetra_param_vec[i] = 
00086   Teuchos::rcp(new Epetra_Vector(*(epetra_param_map[i])));
00087   
00088       // Set parameters
00089       for (unsigned int j=0; j<sacado_param_vec[i]->size(); j++)
00090   (*(epetra_param_vec[i]))[j] = (*(sacado_param_vec[i]))[j].baseValue;
00091 
00092     }
00093 
00094     // Create storage for SG parameter values
00095 #if SG_ACTIVE
00096     if (supports_sg)
00097       p_sg_vals.resize(sg_param_names->size());
00098 #endif
00099   }
00100 
00101   supports_g = (app->getResponseMap() != Teuchos::null);
00102                   
00103 }
00104 
00105 // Overridden from EpetraExt::ModelEvaluator
00106 
00107 Teuchos::RCP<const Epetra_Map>
00108 FEApp::ModelEvaluator::get_x_map() const
00109 {
00110   return app->getMap();
00111 }
00112 
00113 Teuchos::RCP<const Epetra_Map>
00114 FEApp::ModelEvaluator::get_f_map() const
00115 {
00116   return app->getMap();
00117 }
00118 
00119 Teuchos::RCP<const Epetra_Map>
00120 FEApp::ModelEvaluator::get_p_map(int l) const
00121 {
00122   TEST_FOR_EXCEPTION(supports_p == false, 
00123                      Teuchos::Exceptions::InvalidParameter,
00124                      std::endl << 
00125                      "Error!  FEApp::ModelEvaluator::get_p_map():  " <<
00126                      "No parameters have been supplied.  " <<
00127                      "Supplied index l = " << l << std::endl);
00128   TEST_FOR_EXCEPTION(l >= static_cast<int>(epetra_param_map.size()) || l < 0, 
00129          Teuchos::Exceptions::InvalidParameter,
00130                      std::endl << 
00131                      "Error!  FEApp::ModelEvaluator::get_p_map():  " <<
00132                      "Invalid parameter index l = " << l << std::endl);
00133 
00134   return epetra_param_map[l];
00135 }
00136 
00137 Teuchos::RCP<const Epetra_Map>
00138 FEApp::ModelEvaluator::get_g_map(int l) const
00139 {
00140   TEST_FOR_EXCEPTION(supports_g == false, 
00141                      Teuchos::Exceptions::InvalidParameter,
00142                      std::endl << 
00143                      "Error!  FEApp::ModelEvaluator::get_g_map():  " <<
00144                      "No response functions have been supplied.  " <<
00145                      "Supplied index l = " << l << std::endl);
00146   TEST_FOR_EXCEPTION(l != 0, Teuchos::Exceptions::InvalidParameter,
00147                      std::endl << 
00148                      "Error!  FEApp::ModelEvaluator::get_g_map() only " <<
00149                      " supports 1 response vector.  Supplied index l = " << 
00150                      l << std::endl);
00151 
00152   return app->getResponseMap();
00153 }
00154 
00155 Teuchos::RCP<const Teuchos::Array<std::string> >
00156 FEApp::ModelEvaluator::get_p_names(int l) const
00157 {
00158   TEST_FOR_EXCEPTION(supports_p == false, 
00159                      Teuchos::Exceptions::InvalidParameter,
00160                      std::endl << 
00161                      "Error!  FEApp::ModelEvaluator::get_p_names():  " <<
00162                      "No parameters have been supplied.  " <<
00163                      "Supplied index l = " << l << std::endl);
00164   TEST_FOR_EXCEPTION(l >= static_cast<int>(param_names.size()) || l < 0, 
00165          Teuchos::Exceptions::InvalidParameter,
00166                      std::endl << 
00167                      "Error!  FEApp::ModelEvaluator::get_p_names():  " <<
00168                      "Invalid parameter index l = " << l << std::endl);
00169 
00170   return param_names[l];
00171 }
00172 
00173 Teuchos::RCP<const Epetra_Vector>
00174 FEApp::ModelEvaluator::get_x_init() const
00175 {
00176   return app->getInitialSolution();
00177 }
00178 
00179 Teuchos::RCP<const Epetra_Vector>
00180 FEApp::ModelEvaluator::get_p_init(int l) const
00181 {
00182   TEST_FOR_EXCEPTION(supports_p == false, 
00183                      Teuchos::Exceptions::InvalidParameter,
00184                      std::endl << 
00185                      "Error!  FEApp::ModelEvaluator::get_p_init():  " <<
00186                      "No parameters have been supplied.  " <<
00187                      "Supplied index l = " << l << std::endl);
00188   TEST_FOR_EXCEPTION(l >= static_cast<int>(param_names.size()) || l < 0, 
00189          Teuchos::Exceptions::InvalidParameter,
00190                      std::endl << 
00191                      "Error!  FEApp::ModelEvaluator::get_p_init():  " <<
00192                      "Invalid parameter index l = " << l << std::endl);
00193   
00194   return epetra_param_vec[l];
00195 }
00196 
00197 Teuchos::RCP<Epetra_Operator>
00198 FEApp::ModelEvaluator::create_W() const
00199 {
00200   my_W = app->createW();
00201   return my_W;
00202 }
00203 
00204 EpetraExt::ModelEvaluator::InArgs
00205 FEApp::ModelEvaluator::createInArgs() const
00206 {
00207   InArgsSetup inArgs;
00208   inArgs.setModelEvalDescription(this->description());
00209   inArgs.setSupports(IN_ARG_x,true);
00210   inArgs.set_Np(param_names.size());
00211   if (supports_sg) {
00212     inArgs.setSupports(IN_ARG_x_sg,true);
00213     inArgs.set_Np_sg(1); // 1 SG parameter vector
00214   }
00215   else
00216     inArgs.set_Np_sg(0);
00217   if (app->isTransient()) {
00218     inArgs.setSupports(IN_ARG_t,true);
00219     inArgs.setSupports(IN_ARG_x_dot,true);
00220     inArgs.setSupports(IN_ARG_alpha,true);
00221     inArgs.setSupports(IN_ARG_beta,true);
00222     if (supports_sg)
00223       inArgs.setSupports(IN_ARG_x_dot_sg,true);
00224   }
00225   
00226   return inArgs;
00227 }
00228 
00229 EpetraExt::ModelEvaluator::OutArgs
00230 FEApp::ModelEvaluator::createOutArgs() const
00231 {
00232   OutArgsSetup outArgs;
00233   outArgs.setModelEvalDescription(this->description());
00234 
00235   if (supports_p && supports_g) {
00236     outArgs.set_Np_Ng(param_names.size(), 1);
00237     for (unsigned int i=0; i<param_names.size(); i++)
00238       outArgs.setSupports(OUT_ARG_DgDp, 0, i, 
00239         DerivativeSupport(DERIV_MV_BY_COL));
00240   }
00241   else if (supports_p)
00242     outArgs.set_Np_Ng(1, 0);
00243   else if (supports_g)
00244     outArgs.set_Np_Ng(0, 1);
00245   else
00246     outArgs.set_Np_Ng(0, 0);
00247 
00248   if (supports_g) {
00249     outArgs.setSupports(OUT_ARG_DgDx, 0, DerivativeSupport(DERIV_MV_BY_COL));
00250     if (app->isTransient())
00251       outArgs.setSupports(OUT_ARG_DgDx_dot, 0, 
00252                           DerivativeSupport(DERIV_MV_BY_COL));
00253   }
00254   if (supports_p)
00255     for (unsigned int i=0; i<param_names.size(); i++)
00256       outArgs.setSupports(OUT_ARG_DfDp, i, DerivativeSupport(DERIV_MV_BY_COL));
00257 
00258   outArgs.setSupports(OUT_ARG_f,true);
00259   outArgs.setSupports(OUT_ARG_W,true);
00260   outArgs.set_W_properties(
00261     DerivativeProperties(
00262       DERIV_LINEARITY_UNKNOWN ,DERIV_RANK_FULL ,true)
00263     );
00264   if (supports_sg) {
00265     outArgs.setSupports(OUT_ARG_f_sg,true);
00266     outArgs.setSupports(OUT_ARG_W_sg,true);
00267     if (supports_g)
00268       outArgs.set_Np_Ng_sg(0, 1);
00269   }
00270 
00271   return outArgs;
00272 }
00273 
00274 void 
00275 FEApp::ModelEvaluator::evalModel(const InArgs& inArgs, 
00276          const OutArgs& outArgs) const
00277 {
00278   //
00279   // Get the input arguments
00280   //
00281   Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
00282   Teuchos::RCP<const Epetra_Vector> x_dot;
00283   double alpha = 0.0;
00284   double beta = 1.0;
00285   if (app->isTransient()) {
00286     x_dot = inArgs.get_x_dot();
00287     if (x_dot != Teuchos::null) {
00288       alpha = inArgs.get_alpha();
00289       beta = inArgs.get_beta();
00290     }
00291   }
00292   for (int i=0; i<inArgs.Np(); i++) {
00293     Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i);
00294     if (p != Teuchos::null) {
00295       for (unsigned int j=0; j<sacado_param_vec[i]->size(); j++)
00296   (*(sacado_param_vec[i]))[j].baseValue = (*p)[j];
00297     }
00298   }
00299 
00300   //
00301   // Get the output arguments
00302   //
00303   EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out = outArgs.get_f();
00304   Teuchos::RCP<Epetra_Operator> W_out = outArgs.get_W();
00305   
00306   //
00307   // Compute the functions
00308   //
00309   bool f_computed = false;
00310 
00311   // W matrix
00312   if (f_out != Teuchos::null && eval_W_with_f) {
00313     app->computeGlobalJacobian(alpha, beta, x_dot.get(), *x, sacado_param_vec,
00314                                  f_out.get(), *my_W);
00315     f_computed = true;
00316   }
00317   else if (W_out != Teuchos::null && !eval_W_with_f) {
00318     if (f_out.getType() == EVAL_TYPE_EXACT ||
00319         f_out.getType() == EVAL_TYPE_APPROX_DERIV)
00320       app->computeGlobalJacobian(alpha, beta, x_dot.get(), *x, sacado_param_vec,
00321                                  f_out.get(), *W_out);
00322     else
00323       app->computeGlobalPreconditioner(alpha, beta, x_dot.get(), *x, 
00324                                        sacado_param_vec, 
00325                f_out.get(), *W_out);
00326     f_computed = true;
00327   }
00328   
00329   // df/dp
00330   if (supports_p) {
00331     for (int i=0; i<outArgs.Np(); i++) {
00332       Teuchos::RCP<Epetra_MultiVector> dfdp_out = 
00333   outArgs.get_DfDp(i).getMultiVector();
00334       if (dfdp_out != Teuchos::null) {
00335   Teuchos::Array<int> p_indexes = 
00336     outArgs.get_DfDp(i).getDerivativeMultiVector().getParamIndexes();
00337   unsigned int n_params = p_indexes.size();
00338   Teuchos::RCP<ParamVec> p_vec;
00339   if (n_params > 0) {
00340     Teuchos::Array<std::string> p_names(n_params);
00341     for (unsigned int j=0; j<n_params; j++)
00342       p_names[j] = (*(param_names[i]))[p_indexes[j]];
00343     p_vec = Teuchos::rcp(new ParamVec);
00344     app->getParamLib()->fillVector<FEApp::ResidualType>(p_names, *p_vec);
00345     for (unsigned int j=0; j<p_vec->size(); j++)
00346       (*p_vec)[j].baseValue = 
00347         (*(sacado_param_vec[i]))[p_indexes[j]].baseValue;
00348   }
00349   else
00350     p_vec = Teuchos::rcp(sacado_param_vec[0].get(), false);
00351   
00352   app->computeGlobalTangent(0.0, 0.0, false, x_dot.get(), *x, 
00353           sacado_param_vec, p_vec.get(),
00354           NULL, NULL, f_out.get(), NULL, 
00355           dfdp_out.get());
00356 
00357   f_computed = true;
00358       }
00359     }
00360   }
00361 
00362   // f
00363   if(f_out != Teuchos::null && !f_computed) {
00364     app->computeGlobalResidual(x_dot.get(), *x, sacado_param_vec, *f_out);
00365   }
00366 
00367   // Response functions
00368   if (outArgs.Ng() > 0 && supports_g) {
00369     Teuchos::RCP<Epetra_Vector> g_out = outArgs.get_g(0);
00370     Teuchos::RCP<Epetra_MultiVector> dgdx_out = 
00371       outArgs.get_DgDx(0).getMultiVector();
00372     Teuchos::RCP<Epetra_MultiVector> dgdxdot_out;
00373     if (app->isTransient())
00374       dgdxdot_out = outArgs.get_DgDx_dot(0).getMultiVector();
00375     
00376     Teuchos::Array< Teuchos::RCP<ParamVec> > p_vec(outArgs.Np());
00377     Teuchos::Array< Teuchos::RCP<Epetra_MultiVector> > dgdp_out(outArgs.Np());
00378     bool have_dgdp = false;
00379     for (int i=0; i<outArgs.Np(); i++) {
00380       dgdp_out[i] = outArgs.get_DgDp(0,i).getMultiVector();
00381       if (dgdp_out[i] != Teuchos::null)
00382   have_dgdp = true;
00383       Teuchos::Array<int> p_indexes = 
00384   outArgs.get_DgDp(0,i).getDerivativeMultiVector().getParamIndexes();
00385       unsigned int n_params = p_indexes.size();
00386       if (n_params > 0) {
00387   Teuchos::Array<std::string> p_names(n_params);
00388   for (unsigned int j=0; j<n_params; j++)
00389     p_names[i] = (*(param_names[i]))[p_indexes[j]];
00390   p_vec[i] = Teuchos::rcp(new ParamVec);
00391   app->getParamLib()->fillVector<FEApp::ResidualType>(p_names, 
00392                   *(p_vec[i]));
00393   for (unsigned int j=0; j<p_vec[i]->size(); j++)
00394     (*(p_vec[i]))[j].baseValue = (*(sacado_param_vec[i]))[p_indexes[j]].baseValue;
00395       }
00396       else
00397   p_vec[i] = sacado_param_vec[i];
00398     }
00399 
00400     if (have_dgdp ||dgdx_out != Teuchos::null || dgdxdot_out != Teuchos::null) {
00401       app->evaluateResponseGradients(x_dot.get(), *x, sacado_param_vec, p_vec, 
00402                                      g_out.get(), dgdx_out.get(), 
00403                                      dgdxdot_out.get(), dgdp_out);
00404     }
00405     else if (g_out != Teuchos::null)
00406       app->evaluateResponses(x_dot.get(), *x, sacado_param_vec, *g_out);
00407   }
00408 
00409 #if SG_ACTIVE
00410   if (supports_sg) {
00411     InArgs::sg_const_vector_t x_sg = 
00412       inArgs.get_x_sg();
00413     if (x_sg != Teuchos::null) {
00414       InArgs::sg_const_vector_t x_dot_sg;
00415       if (app->isTransient())
00416   x_dot_sg = inArgs.get_x_dot_sg();
00417       InArgs::sg_const_vector_t epetra_p_sg = 
00418   inArgs.get_p_sg(0);
00419       Teuchos::Array<SGType> *p_sg_ptr = NULL;
00420       if (epetra_p_sg != Teuchos::null) {
00421   for (unsigned int i=0; i<p_sg_vals.size(); i++) {
00422     int num_sg_blocks = epetra_p_sg->size();
00423     p_sg_vals[i].reset(app->getStochasticExpansion());
00424     p_sg_vals[i].copyForWrite();
00425     for (int j=0; j<num_sg_blocks; j++) {
00426       p_sg_vals[i].fastAccessCoeff(j) = (*epetra_p_sg)[j][i];
00427     }
00428   }
00429   p_sg_ptr = &p_sg_vals;
00430       }
00431       OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
00432       OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
00433       if (W_sg != Teuchos::null)
00434   app->computeGlobalSGJacobian(alpha, beta, x_dot_sg.get(), *x_sg, 
00435              sacado_param_vec[0].get(), 
00436              sacado_param_vec[1].get(), p_sg_ptr,
00437              f_sg.get(), *W_sg);
00438       else if (f_sg != Teuchos::null)
00439   app->computeGlobalSGResidual(x_dot_sg.get(), *x_sg, 
00440              sacado_param_vec[0].get(), 
00441              sacado_param_vec[1].get(), p_sg_ptr,
00442              *f_sg);
00443 
00444       // Response functions
00445       if (outArgs.Ng() > 0 && supports_g) {
00446   OutArgs::sg_vector_t g_sg = outArgs.get_g_sg(0);
00447   if (g_sg != Teuchos::null)
00448     app->evaluateSGResponses(x_dot_sg.get(), *x_sg, 
00449            sacado_param_vec[0].get(), 
00450            sacado_param_vec[1].get(), p_sg_ptr, 
00451            *g_sg);
00452       }
00453     }
00454   }
00455 #endif
00456 }

Generated on Wed May 12 21:39:33 2010 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.4.7