FEApp_Application.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_Application.hpp"
00033 #include "FEApp_ProblemFactory.hpp"
00034 #include "FEApp_QuadratureFactory.hpp"
00035 #include "FEApp_DiscretizationFactory.hpp"
00036 #include "Epetra_LocalMap.h"
00037 #if SG_ACTIVE
00038 #include "FEApp_SGGaussQuadResidualGlobalFill.hpp"
00039 #include "FEApp_SGGaussQuadJacobianGlobalFill.hpp"
00040 #include "Stokhos_MatrixFreeEpetraOp.hpp"
00041 #include "Stokhos_MeanEpetraOp.hpp"
00042 #endif
00043 #include "Teuchos_TimeMonitor.hpp"
00044 
00045 FEApp::Application::Application(
00046        const std::vector<double>& coords,
00047        const Teuchos::RCP<const Epetra_Comm>& comm,
00048        const Teuchos::RCP<Teuchos::ParameterList>& params_,
00049        bool is_transient,
00050        const Epetra_Vector* initial_soln) :
00051   params(params_),
00052   transient(is_transient)
00053 {
00054   // Create parameter library
00055   paramLib = Teuchos::rcp(new ParamLib);
00056 
00057   // Create problem object
00058   Teuchos::RCP<Teuchos::ParameterList> problemParams = 
00059     Teuchos::rcp(&(params->sublist("Problem")),false);
00060   FEApp::ProblemFactory problemFactory(problemParams, paramLib);
00061   Teuchos::RCP<FEApp::AbstractProblem> problem = 
00062     problemFactory.create();
00063 
00064   // Get number of equations
00065   unsigned int num_equations = problem->numEquations();
00066 
00067   // Create quadrature object
00068   Teuchos::RCP<Teuchos::ParameterList> quadParams = 
00069     Teuchos::rcp(&(params->sublist("Quadrature")),false);
00070   FEApp::QuadratureFactory quadFactory(quadParams);
00071   quad = quadFactory.create();
00072 
00073   // Create discretization object
00074   Teuchos::RCP<Teuchos::ParameterList> discParams = 
00075     Teuchos::rcp(&(params->sublist("Discretization")),false);
00076   FEApp::DiscretizationFactory discFactory(discParams);
00077   disc = discFactory.create(coords, num_equations, comm);
00078   disc->createMesh();
00079   disc->createMaps();
00080   disc->createJacobianGraphs();
00081 
00082   // Create Epetra objects
00083   importer = Teuchos::rcp(new Epetra_Import(*(disc->getOverlapMap()), 
00084                                             *(disc->getMap())));
00085   exporter = Teuchos::rcp(new Epetra_Export(*(disc->getOverlapMap()), 
00086                                             *(disc->getMap())));
00087   overlapped_x = Teuchos::rcp(new Epetra_Vector(*(disc->getOverlapMap())));
00088   if (transient)
00089     overlapped_xdot = 
00090       Teuchos::rcp(new Epetra_Vector(*(disc->getOverlapMap())));
00091   overlapped_f = Teuchos::rcp(new Epetra_Vector(*(disc->getOverlapMap())));
00092   overlapped_jac = 
00093     Teuchos::rcp(new Epetra_CrsMatrix(Copy, 
00094                                       *(disc->getOverlapJacobianGraph())));
00095 
00096   // Initialize problem
00097   initial_x = Teuchos::rcp(new Epetra_Vector(*(disc->getMap())));
00098   problem->buildProblem(*(disc->getMap()), *(disc->getOverlapMap()), 
00099                         pdeTM, bc, responses, initial_x);
00100   if (initial_soln != NULL)
00101     initial_x->Scale(1.0, *initial_soln);
00102   typedef FEApp::AbstractPDE_TemplateManager<EvalTypes>::iterator iterator;
00103   int nqp = quad->numPoints();
00104   int nn = disc->getNumNodesPerElement();
00105   for (iterator it = pdeTM.begin(); it != pdeTM.end(); ++it)
00106     it->init(nqp, nn);
00107 
00108   // Create response map
00109   unsigned int total_num_responses = 0;
00110   for (unsigned int i=0; i<responses.size(); i++)
00111     total_num_responses += responses[i]->numResponses();
00112   if (total_num_responses > 0)
00113     response_map = Teuchos::rcp(new Epetra_LocalMap(total_num_responses, 0,
00114                                                     *comm));
00115 
00116 #if SG_ACTIVE
00117   bool enable_sg = params->get("Enable Stochastic Galerkin",false);
00118   if (enable_sg) {
00119     sg_expansion = params->get< Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > >("Stochastic Galerkin expansion");
00120     sg_quad = params->get< Teuchos::RCP<const Stokhos::Quadrature<int,double> > >("Stochastic Galerkin quadrature");
00121     sg_basis = sg_expansion->getBasis();
00122 
00123     // Create Epetra orthogonal polynomial objects
00124     sg_overlapped_x = 
00125       Teuchos::rcp(new Stokhos::VectorOrthogPoly<Epetra_Vector>(sg_basis,
00126                 *overlapped_x));
00127     if (transient)
00128       sg_overlapped_xdot = 
00129   Teuchos::rcp(new Stokhos::VectorOrthogPoly<Epetra_Vector>(sg_basis,
00130                   *overlapped_xdot));
00131     sg_overlapped_f = 
00132       Teuchos::rcp(new Stokhos::VectorOrthogPoly<Epetra_Vector>(sg_basis,
00133                 *overlapped_f));
00134     sg_overlapped_jac = 
00135       Teuchos::rcp(new Stokhos::VectorOrthogPoly<Epetra_CrsMatrix>(sg_basis,
00136                    *overlapped_jac));
00137   }
00138 #endif
00139 }
00140 
00141 FEApp::Application::~Application()
00142 {
00143 }
00144 
00145 Teuchos::RCP<const Epetra_Map>
00146 FEApp::Application::getMap() const
00147 {
00148   return disc->getMap();
00149 }
00150 
00151 Teuchos::RCP<const Epetra_Map>
00152 FEApp::Application::getResponseMap() const
00153 {
00154   return response_map;
00155 }
00156 
00157 Teuchos::RCP<const Epetra_CrsGraph>
00158 FEApp::Application::getJacobianGraph() const
00159 {
00160   return disc->getJacobianGraph();
00161 }
00162 
00163 Teuchos::RCP<const Epetra_Vector>
00164 FEApp::Application::getInitialSolution() const
00165 {
00166   return initial_x;
00167 }
00168 
00169 Teuchos::RCP<ParamLib> 
00170 FEApp::Application::getParamLib()
00171 {
00172   return paramLib;
00173 }
00174 
00175 bool
00176 FEApp::Application::isTransient() const
00177 {
00178   return transient;
00179 }
00180 
00181 Teuchos::RCP<Epetra_Operator>
00182 FEApp::Application::createW() const
00183 {
00184   return Teuchos::rcp(new Epetra_CrsMatrix(::Copy, 
00185              *(disc->getJacobianGraph())));
00186 }
00187 
00188 Teuchos::RCP<Epetra_Operator>
00189 FEApp::Application::createPrec() const
00190 {
00191   return Teuchos::rcp(new Epetra_CrsMatrix(::Copy, 
00192              *(disc->getJacobianGraph())));
00193 }
00194 
00195 void
00196 FEApp::Application::computeGlobalResidual(
00197                           const Epetra_Vector* xdot,
00198         const Epetra_Vector& x,
00199         const Teuchos::Array< Teuchos::RCP<ParamVec> >& p,
00200         Epetra_Vector& f)
00201 {
00202   TEUCHOS_FUNC_TIME_MONITOR("FEApp::Application::computeGlobalResidual");
00203 
00204   // Scatter x to the overlapped distrbution
00205   overlapped_x->Import(x, *importer, Insert);
00206 
00207   // Scatter xdot to the overlapped distribution
00208   if (transient)
00209     overlapped_xdot->Import(*xdot, *importer, Insert);
00210 
00211   // Set parameters
00212   for (unsigned int i=0; i<p.size(); i++) {
00213     if (p[i] != Teuchos::null)
00214       for (unsigned int j=0; j<p[i]->size(); j++)
00215   (*(p[i]))[j].family->setRealValueForAllTypes((*(p[i]))[j].baseValue);
00216   }
00217 
00218   // Zero out overlapped residual
00219   overlapped_f->PutScalar(0.0);
00220   f.PutScalar(0.0);
00221 
00222   // Create residual init/post op
00223   Teuchos::RCP<FEApp::ResidualOp> op = 
00224     Teuchos::rcp(new FEApp::ResidualOp(overlapped_xdot, overlapped_x, 
00225                                        overlapped_f));
00226 
00227   // Get template PDE instantiation
00228   Teuchos::RCP< FEApp::AbstractPDE<FEApp::ResidualType> > pde = 
00229     pdeTM.getAsObject<FEApp::ResidualType>();
00230 
00231   // Do global fill
00232   FEApp::GlobalFill<FEApp::ResidualType> globalFill(disc->getMesh(), quad, 
00233                                                     pde, bc, transient);
00234   globalFill.computeGlobalFill(*op);
00235 
00236   // Assemble global residual
00237   f.Export(*overlapped_f, *exporter, Add);
00238 }
00239 
00240 void
00241 FEApp::Application::computeGlobalJacobian(
00242           double alpha, double beta,
00243           const Epetra_Vector* xdot,
00244           const Epetra_Vector& x,
00245           const Teuchos::Array< Teuchos::RCP<ParamVec> >& p,
00246           Epetra_Vector* f,
00247           Epetra_Operator& jacOp)
00248 {
00249   TEUCHOS_FUNC_TIME_MONITOR("FEApp::Application::computeGlobalJacobian");
00250 
00251   // Cast jacOp to an Epetra_CrsMatrix
00252   Epetra_CrsMatrix& jac = dynamic_cast<Epetra_CrsMatrix&>(jacOp);
00253 
00254   // Scatter x to the overlapped distrbution
00255   overlapped_x->Import(x, *importer, Insert);
00256 
00257   // Scatter xdot to the overlapped distribution
00258   if (transient)
00259     overlapped_xdot->Import(*xdot, *importer, Insert);
00260 
00261   // Set parameters
00262   for (unsigned int i=0; i<p.size(); i++) {
00263     if (p[i] != Teuchos::null)
00264       for (unsigned int j=0; j<p[i]->size(); j++)
00265   (*(p[i]))[j].family->setRealValueForAllTypes((*(p[i]))[j].baseValue);
00266   }
00267 
00268   // Zero out overlapped residual
00269   Teuchos::RCP<Epetra_Vector> overlapped_ff;
00270   if (f != NULL) {
00271     overlapped_ff = overlapped_f;
00272     overlapped_ff->PutScalar(0.0);
00273     f->PutScalar(0.0);
00274   }
00275 
00276   // Zero out Jacobian
00277   overlapped_jac->PutScalar(0.0);
00278   jac.PutScalar(0.0);
00279 
00280   // Create Jacobian init/post op
00281   Teuchos::RCP<FEApp::JacobianOp> op
00282     = Teuchos::rcp(new FEApp::JacobianOp(alpha, beta, overlapped_xdot, 
00283                                          overlapped_x, overlapped_ff, 
00284                                          overlapped_jac));
00285 
00286   // Get template PDE instantiation
00287   Teuchos::RCP< FEApp::AbstractPDE<FEApp::JacobianType> > pde = 
00288     pdeTM.getAsObject<FEApp::JacobianType>();
00289 
00290   // Do global fill
00291   FEApp::GlobalFill<FEApp::JacobianType> globalFill(disc->getMesh(), 
00292                                                     quad, pde, bc, 
00293                                                     transient);
00294   globalFill.computeGlobalFill(*op);
00295 
00296   // Assemble global residual
00297   if (f != NULL)
00298     f->Export(*overlapped_f, *exporter, Add);
00299 
00300   // Assemble global Jacobian
00301   jac.Export(*overlapped_jac, *exporter, Add);
00302 
00303   jac.FillComplete(true);
00304 }
00305 
00306 void
00307 FEApp::Application::computeGlobalPreconditioner(
00308           double alpha, double beta,
00309           const Epetra_Vector* xdot,
00310           const Epetra_Vector& x,
00311           const Teuchos::Array< Teuchos::RCP<ParamVec> >& p,
00312           Epetra_Vector* f,
00313           Epetra_Operator& jacOp)
00314 {
00315   TEUCHOS_FUNC_TIME_MONITOR("FEApp::Application::computeGlobalPreconditioner");
00316 
00317   // Cast jacOp to an Epetra_CrsMatrix
00318   Epetra_CrsMatrix& jac = dynamic_cast<Epetra_CrsMatrix&>(jacOp);
00319 
00320   // Scatter x to the overlapped distrbution
00321   overlapped_x->Import(x, *importer, Insert);
00322 
00323   // Scatter xdot to the overlapped distribution
00324   if (transient)
00325     overlapped_xdot->Import(*xdot, *importer, Insert);
00326 
00327   // Set parameters
00328   for (unsigned int i=0; i<p.size(); i++) {
00329     if (p[i] != Teuchos::null)
00330       for (unsigned int j=0; j<p[i]->size(); j++)
00331   (*(p[i]))[j].family->setRealValueForAllTypes((*(p[i]))[j].baseValue);
00332   }
00333 
00334   // Zero out overlapped residual
00335   Teuchos::RCP<Epetra_Vector> overlapped_ff;
00336   if (f != NULL) {
00337     overlapped_ff = overlapped_f;
00338     overlapped_ff->PutScalar(0.0);
00339     f->PutScalar(0.0);
00340   }
00341 
00342   // Zero out Jacobian
00343   overlapped_jac->PutScalar(0.0);
00344   jac.PutScalar(0.0);
00345 
00346   // Create Jacobian init/post op
00347   Teuchos::RCP<FEApp::JacobianOp> op = 
00348     Teuchos::rcp(new FEApp::JacobianOp(alpha, beta, overlapped_xdot, 
00349                                        overlapped_x, overlapped_ff, 
00350                                        overlapped_jac));
00351 
00352   // Get template PDE instantiation
00353   Teuchos::RCP< FEApp::AbstractPDE<FEApp::JacobianType> > pde = 
00354     pdeTM.getAsObject<FEApp::JacobianType>();
00355 
00356   // Do global fill
00357   FEApp::GlobalFill<FEApp::JacobianType> globalFill(disc->getMesh(), 
00358                                                     quad, pde, bc, 
00359                                                     transient);
00360   globalFill.computeGlobalFill(*op);
00361 
00362   // Assemble global residual
00363   if (f != NULL)
00364     f->Export(*overlapped_f, *exporter, Add);
00365 
00366   // Assemble global Jacobian
00367   jac.Export(*overlapped_jac, *exporter, Add);
00368 
00369   jac.FillComplete(true);
00370 }
00371 
00372 void
00373 FEApp::Application::computeGlobalTangent(
00374             double alpha, double beta,
00375             bool sum_derivs,
00376             const Epetra_Vector* xdot,
00377             const Epetra_Vector& x,
00378             const Teuchos::Array< Teuchos::RCP<ParamVec> >& p,
00379             ParamVec* deriv_p,
00380             const Epetra_MultiVector* Vx,
00381             const Teuchos::SerialDenseMatrix<int,double>* Vp,
00382             Epetra_Vector* f,
00383             Epetra_MultiVector* JVx,
00384             Epetra_MultiVector* fVp)
00385 {
00386   // Scatter x to the overlapped distrbution
00387   overlapped_x->Import(x, *importer, Insert);
00388 
00389   // Scatter xdot to the overlapped distribution
00390   if (transient)
00391     overlapped_xdot->Import(*xdot, *importer, Insert);
00392 
00393   // Set parameters
00394   for (unsigned int i=0; i<p.size(); i++) {
00395     if (p[i] != Teuchos::null)
00396       for (unsigned int j=0; j<p[i]->size(); j++)
00397   (*(p[i]))[j].family->setRealValueForAllTypes((*(p[i]))[j].baseValue);
00398   }
00399 
00400   // Zero out overlapped residual
00401   Teuchos::RCP<Epetra_Vector> overlapped_ff;
00402   if (f != NULL) {
00403     overlapped_ff = overlapped_f;
00404     overlapped_ff->PutScalar(0.0);
00405     f->PutScalar(0.0);
00406   }
00407 
00408   Teuchos::RCP<Epetra_MultiVector> overlapped_JVx;
00409   if (JVx != NULL) {
00410     overlapped_JVx = 
00411       Teuchos::rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), 
00412                                           JVx->NumVectors()));
00413     overlapped_JVx->PutScalar(0.0);
00414     JVx->PutScalar(0.0);
00415   }
00416   
00417   Teuchos::RCP<Epetra_MultiVector> overlapped_fVp;
00418   if (fVp != NULL) {
00419     overlapped_fVp = 
00420       Teuchos::rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), 
00421                                           fVp->NumVectors()));
00422     overlapped_fVp->PutScalar(0.0);
00423     fVp->PutScalar(0.0);
00424   }
00425 
00426   Teuchos::RCP<Epetra_MultiVector> overlapped_Vx;
00427   if (Vx != NULL) {
00428     overlapped_Vx = 
00429       Teuchos::rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), 
00430                                           Vx->NumVectors()));
00431   }
00432 
00433   Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,double> > vp =
00434     Teuchos::rcp(Vp, false);
00435   Teuchos::RCP<ParamVec> params = 
00436     Teuchos::rcp(deriv_p, false);
00437 
00438   // Create Jacobian init/post op
00439   Teuchos::RCP<FEApp::TangentOp> op = 
00440     Teuchos::rcp(new FEApp::TangentOp(alpha, beta, sum_derivs,
00441                                       overlapped_xdot, 
00442                                       overlapped_x,
00443                                       params,
00444                                       overlapped_Vx,
00445                                       overlapped_Vx,
00446                                       vp,
00447                                       overlapped_ff, 
00448                                       overlapped_JVx,
00449                                       overlapped_fVp));
00450 
00451   // Get template PDE instantiation
00452   Teuchos::RCP< FEApp::AbstractPDE<FEApp::TangentType> > pde = 
00453     pdeTM.getAsObject<FEApp::TangentType>();
00454 
00455   // Do global fill
00456   FEApp::GlobalFill<FEApp::TangentType> globalFill(disc->getMesh(), 
00457                                                    quad, pde, bc, 
00458                                                    transient);
00459   globalFill.computeGlobalFill(*op);
00460 
00461   // Assemble global residual
00462   if (f != NULL)
00463     f->Export(*overlapped_f, *exporter, Add);
00464 
00465   // Assemble derivatives
00466   if (JVx != NULL)
00467     JVx->Export(*overlapped_JVx, *exporter, Add);
00468   if (fVp != NULL)
00469     fVp->Export(*overlapped_fVp, *exporter, Add);
00470 }
00471 
00472 void
00473 FEApp::Application::
00474 evaluateResponses(const Epetra_Vector* xdot,
00475                   const Epetra_Vector& x,
00476                   const Teuchos::Array< Teuchos::RCP<ParamVec> >& p,
00477                   Epetra_Vector& g)
00478 {
00479   const Epetra_Comm& comm = x.Map().Comm();
00480   unsigned int offset = 0;
00481   for (unsigned int i=0; i<responses.size(); i++) {
00482 
00483     // Create Epetra_Map for response function
00484     unsigned int num_responses = responses[i]->numResponses();
00485     Epetra_LocalMap local_response_map(num_responses, 0, comm);
00486 
00487     // Create Epetra_Vector for response function
00488     Epetra_Vector local_g(local_response_map);
00489 
00490     // Evaluate response function
00491     responses[i]->evaluateResponses(xdot, x, p, local_g);
00492 
00493     // Copy result into combined result
00494     for (unsigned int j=0; j<num_responses; j++)
00495       g[offset+j] = local_g[j];
00496 
00497     // Increment offset in combined result
00498     offset += num_responses;
00499   }
00500 }
00501 
00502 void
00503 FEApp::Application::
00504 evaluateResponseTangents(
00505      const Epetra_Vector* xdot,
00506      const Epetra_Vector& x,
00507      const Teuchos::Array< Teuchos::RCP<ParamVec> >& p,
00508      const Teuchos::Array< Teuchos::RCP<ParamVec> >& deriv_p,
00509      const Teuchos::Array< Teuchos::RCP<Epetra_MultiVector> >& dxdot_dp,
00510      const Teuchos::Array< Teuchos::RCP<Epetra_MultiVector> >& dx_dp,
00511      Epetra_Vector* g,
00512      const Teuchos::Array< Teuchos::RCP<Epetra_MultiVector> >& gt)
00513 {
00514   const Epetra_Comm& comm = x.Map().Comm();
00515   unsigned int offset = 0;
00516   Teuchos::Array< Teuchos::RCP<Epetra_MultiVector> > local_gt(gt.size());
00517   for (unsigned int i=0; i<responses.size(); i++) {
00518 
00519     // Create Epetra_Map for response function
00520     unsigned int num_responses = responses[i]->numResponses();
00521     Epetra_LocalMap local_response_map(num_responses, 0, comm);
00522 
00523     // Create Epetra_Vectors for response function
00524     Teuchos::RCP<Epetra_Vector> local_g;
00525     if (g != NULL)
00526       local_g = Teuchos::rcp(new Epetra_Vector(local_response_map));
00527     for (unsigned int j=0; j<gt.size(); j++)
00528       if (gt[j] != Teuchos::null)
00529   local_gt[j] = Teuchos::rcp(new Epetra_MultiVector(local_response_map, 
00530                 gt[j]->NumVectors()));
00531 
00532     // Evaluate response function
00533     responses[i]->evaluateTangents(xdot, x, p, deriv_p, dxdot_dp, dx_dp, 
00534            local_g.get(), local_gt);
00535 
00536     // Copy results into combined result
00537     for (unsigned int j=0; j<num_responses; j++) {
00538       if (g != NULL)
00539         (*g)[offset+j] = (*local_g)[j];
00540       for (unsigned int l=0; l<gt.size(); l++)
00541   if (gt[l] != Teuchos::null)
00542     for (int k=0; k<gt[l]->NumVectors(); k++)
00543       (*gt[l])[k][offset+j] = (*local_gt[l])[k][j];
00544     }
00545 
00546     // Increment offset in combined result
00547     offset += num_responses;
00548   }
00549 }
00550 
00551 void
00552 FEApp::Application::
00553 evaluateResponseGradients(
00554       const Epetra_Vector* xdot,
00555       const Epetra_Vector& x,
00556       const Teuchos::Array< Teuchos::RCP<ParamVec> >& p,
00557       const Teuchos::Array< Teuchos::RCP<ParamVec> >& deriv_p,
00558       Epetra_Vector* g,
00559       Epetra_MultiVector* dg_dx,
00560       Epetra_MultiVector* dg_dxdot,
00561       const Teuchos::Array< Teuchos::RCP<Epetra_MultiVector> >& dg_dp)
00562 {
00563   const Epetra_Comm& comm = x.Map().Comm();
00564   unsigned int offset = 0;
00565   Teuchos::Array< Teuchos::RCP<Epetra_MultiVector> > local_dgdp(dg_dp.size());
00566   for (unsigned int i=0; i<responses.size(); i++) {
00567 
00568     // Create Epetra_Map for response function
00569     unsigned int num_responses = responses[i]->numResponses();
00570     Epetra_LocalMap local_response_map(num_responses, 0, comm);
00571 
00572     // Create Epetra_Vectors for response function
00573     Teuchos::RCP<Epetra_Vector> local_g;
00574     if (g != NULL)
00575       local_g = Teuchos::rcp(new Epetra_Vector(local_response_map));
00576     Teuchos::RCP<Epetra_MultiVector> local_dgdx;
00577     if (dg_dx != NULL)
00578       local_dgdx = Teuchos::rcp(new Epetra_MultiVector(dg_dx->Map(), 
00579                                                        num_responses));
00580     Teuchos::RCP<Epetra_MultiVector> local_dgdxdot;
00581     if (dg_dxdot != NULL)
00582       local_dgdxdot = Teuchos::rcp(new Epetra_MultiVector(dg_dxdot->Map(), 
00583                                                           num_responses));
00584 
00585     for (unsigned int j=0; j<dg_dp.size(); j++)
00586       if (dg_dp[j] != Teuchos::null)
00587   local_dgdp[j] = Teuchos::rcp(new Epetra_MultiVector(local_response_map, 
00588                   dg_dp[j]->NumVectors()));
00589 
00590     // Evaluate response function
00591     responses[i]->evaluateGradients(xdot, x, p, deriv_p, local_g.get(), 
00592                                     local_dgdx.get(), local_dgdxdot.get(), 
00593                                     local_dgdp);
00594 
00595     // Copy results into combined result
00596     for (unsigned int j=0; j<num_responses; j++) {
00597       if (g != NULL)
00598         (*g)[offset+j] = (*local_g)[j];
00599       if (dg_dx != NULL)
00600         (*dg_dx)(offset+j)->Update(1.0, *((*local_dgdx)(j)), 0.0);
00601       if (dg_dxdot != NULL)
00602         (*dg_dxdot)(offset+j)->Update(1.0, *((*local_dgdxdot)(j)), 0.0);
00603       for (unsigned int l=0; l<dg_dp.size(); l++)
00604   if (dg_dp[l] != Teuchos::null)
00605     for (int k=0; k<dg_dp[l]->NumVectors(); k++)
00606       (*dg_dp[l])[k][offset+j] = (*local_dgdp[l])[k][j];
00607     }
00608 
00609     // Increment offset in combined result
00610     offset += num_responses;
00611   }
00612 }
00613 
00614 #if SG_ACTIVE
00615 void
00616 FEApp::Application::computeGlobalSGResidual(
00617       const Stokhos::VectorOrthogPoly<Epetra_Vector>* sg_xdot,
00618       const Stokhos::VectorOrthogPoly<Epetra_Vector>& sg_x,
00619       const ParamVec* p,
00620       const ParamVec* sg_p,
00621       const Teuchos::Array<SGType>* sg_p_vals,
00622       Stokhos::VectorOrthogPoly<Epetra_Vector>& sg_f)
00623 {
00624   TEUCHOS_FUNC_TIME_MONITOR("FEApp::Application::computeGlobalSGResidual");
00625 
00626   for (int i=0; i<sg_x.size(); i++) {
00627 
00628     // Scatter x to the overlapped distrbution
00629     (*sg_overlapped_x)[i].Import(sg_x[i], *importer, Insert);
00630 
00631     // Scatter xdot to the overlapped distribution
00632     if (transient)
00633       (*sg_overlapped_xdot)[i].Import((*sg_xdot)[i], *importer, Insert);
00634 
00635     // Zero out overlapped residual
00636     (*sg_overlapped_f)[i].PutScalar(0.0);
00637     sg_f[i].PutScalar(0.0);
00638 
00639   }
00640 
00641   // Set real parameters
00642   if (p != NULL) {
00643     for (unsigned int i=0; i<p->size(); ++i) {
00644       (*p)[i].family->setRealValueForAllTypes((*p)[i].baseValue);
00645     }
00646   }
00647 
00648   // Set SG parameters
00649   if (sg_p != NULL && sg_p_vals != NULL) {
00650     for (unsigned int i=0; i<sg_p->size(); ++i) {
00651       (*sg_p)[i].family->setValue<FEApp::SGResidualType>((*sg_p_vals)[i]);
00652     }
00653   }
00654 
00655   // Create residual init/post op
00656   Teuchos::RCP<FEApp::SGResidualOp> sg_res_fill_op = 
00657     Teuchos::rcp(new FEApp::SGResidualOp(sg_expansion, 
00658            sg_overlapped_xdot, 
00659            sg_overlapped_x, 
00660            sg_overlapped_f));
00661     
00662   // Get template PDE instantiation
00663   Teuchos::RCP< FEApp::AbstractPDE<FEApp::SGResidualType> > pde = 
00664     pdeTM.getAsObject<FEApp::SGResidualType>();
00665 
00666   // Instantiate global fill
00667   if (sg_res_global_fill == Teuchos::null) {
00668     std::string method = params->get("SG Method", "AD");
00669     if (method == "AD") {
00670       sg_res_global_fill = 
00671         Teuchos::rcp(new FEApp::GlobalFill<FEApp::SGResidualType>(
00672     disc->getMesh(), quad, pde, bc, transient));
00673     }
00674     else if (method == "Gauss Quadrature") {
00675       Teuchos::RCP< FEApp::AbstractPDE<FEApp::ResidualType> > res_pde = 
00676   pdeTM.getAsObject<FEApp::ResidualType>();
00677       sg_res_global_fill = 
00678         Teuchos::rcp(new FEApp::SGGaussQuadResidualGlobalFill(disc->getMesh(), 
00679                                                               quad, pde, bc, 
00680                                                               transient,
00681                                                               sg_basis,
00682                     sg_quad,
00683                                                               res_pde,
00684                                                               sg_p));
00685     }
00686   }
00687 
00688   // Do global fill
00689   sg_res_global_fill->computeGlobalFill(*sg_res_fill_op);
00690 
00691   // Assemble global residual
00692   for (int i=0; i<sg_f.size(); i++)
00693     sg_f[i].Export((*sg_overlapped_f)[i], *exporter, Add);
00694 }
00695 
00696 void
00697 FEApp::Application::computeGlobalSGJacobian(
00698       double alpha, double beta,
00699       const Stokhos::VectorOrthogPoly<Epetra_Vector>* sg_xdot,
00700       const Stokhos::VectorOrthogPoly<Epetra_Vector>& sg_x,
00701       const ParamVec* p,
00702       const ParamVec* sg_p,
00703       const Teuchos::Array<SGType>* sg_p_vals,
00704       Stokhos::VectorOrthogPoly<Epetra_Vector>* sg_f,
00705       Stokhos::VectorOrthogPoly<Epetra_Operator>& sg_jac)
00706 {
00707   TEUCHOS_FUNC_TIME_MONITOR("FEApp::Application::computeGlobalSGJacobian");
00708 
00709   for (int i=0; i<sg_x.size(); i++) {
00710 
00711     // Scatter x to the overlapped distrbution
00712     (*sg_overlapped_x)[i].Import(sg_x[i], *importer, Insert);
00713 
00714     // Scatter xdot to the overlapped distribution
00715     if (transient)
00716       (*sg_overlapped_xdot)[i].Import((*sg_xdot)[i], *importer, Insert);
00717 
00718     // Zero out overlapped residual
00719     if (sg_f != NULL) {
00720       (*sg_overlapped_f)[i].PutScalar(0.0);
00721       (*sg_f)[i].PutScalar(0.0);
00722     }
00723     (*sg_overlapped_jac)[i].PutScalar(0.0);
00724 
00725   }
00726 
00727   // Set real parameters
00728   if (p != NULL) {
00729     for (unsigned int i=0; i<p->size(); ++i) {
00730       (*p)[i].family->setRealValueForAllTypes((*p)[i].baseValue);
00731     }
00732   }
00733 
00734   // Set SG parameters
00735   if (sg_p != NULL && sg_p_vals != NULL) {
00736     for (unsigned int i=0; i<sg_p->size(); ++i) {
00737       (*sg_p)[i].family->setValue<FEApp::SGJacobianType>((*sg_p_vals)[i]);
00738     }
00739   }
00740 
00741   // Create Jacobian init/post op
00742   Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_Vector> > sg_overlapped_ff;
00743   if (sg_f != NULL)
00744     sg_overlapped_ff = sg_overlapped_f;
00745   Teuchos::RCP<FEApp::SGJacobianOp> sg_jac_fill_op = 
00746     Teuchos::rcp(new FEApp::SGJacobianOp(sg_expansion,
00747            alpha, beta, 
00748            sg_overlapped_xdot, 
00749            sg_overlapped_x, 
00750            sg_overlapped_ff, 
00751            sg_overlapped_jac));
00752 
00753   // Get template PDE instantiation
00754   Teuchos::RCP< FEApp::AbstractPDE<FEApp::SGJacobianType> > pde = 
00755     pdeTM.getAsObject<FEApp::SGJacobianType>();
00756 
00757   // Instantiate global fill
00758   if (sg_jac_global_fill == Teuchos::null) {
00759     std::string method = params->get("SG Method", "AD");
00760     if (method == "AD") {
00761       sg_jac_global_fill = 
00762   Teuchos::rcp(new FEApp::GlobalFill<FEApp::SGJacobianType>(
00763       disc->getMesh(), quad, pde, bc,transient));
00764     }
00765     else if (method == "Gauss Quadrature") {
00766       Teuchos::RCP< FEApp::AbstractPDE<FEApp::JacobianType> > jac_pde = 
00767   pdeTM.getAsObject<FEApp::JacobianType>();
00768       sg_jac_global_fill = 
00769   Teuchos::rcp(new FEApp::SGGaussQuadJacobianGlobalFill(disc->getMesh(),
00770                     quad, pde, bc, 
00771                     transient,
00772                     sg_basis,
00773                     sg_quad,
00774                     jac_pde,
00775                     sg_p, 
00776                     alpha, beta));
00777     }
00778   }
00779 
00780   // Do global fill
00781   sg_jac_global_fill->computeGlobalFill(*sg_jac_fill_op);
00782   
00783   // Assemble global residual
00784   if (sg_f != NULL)
00785     for (int i=0; i<sg_f->size(); i++)
00786       (*sg_f)[i].Export((*sg_overlapped_f)[i], *exporter, Add);
00787     
00788   // Assemble block Jacobians
00789   Teuchos::RCP<Epetra_CrsMatrix> jac;
00790   for (int i=0; i<sg_jac.size(); i++) {
00791     jac = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_jac.getCoeffPtr(i), 
00792                   true);
00793     jac->PutScalar(0.0);
00794     jac->Export((*sg_overlapped_jac)[i], *exporter, Add);
00795     jac->FillComplete(true);
00796   }
00797 }
00798 
00799 void
00800 FEApp::Application::
00801 evaluateSGResponses(const Stokhos::VectorOrthogPoly<Epetra_Vector>* sg_xdot,
00802         const Stokhos::VectorOrthogPoly<Epetra_Vector>& sg_x,
00803         const ParamVec* p,
00804         const ParamVec* sg_p,
00805         const Teuchos::Array<SGType>* sg_p_vals,
00806         Stokhos::VectorOrthogPoly<Epetra_Vector>& sg_g)
00807 {
00808   TEUCHOS_FUNC_TIME_MONITOR("FEApp::Application::evaluateSGResponses");
00809 
00810   const Epetra_Comm& comm = sg_x[0].Map().Comm();
00811   unsigned int offset = 0;
00812   Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis = 
00813     sg_x.basis();
00814   for (unsigned int i=0; i<responses.size(); i++) {
00815 
00816     // Create Epetra_Map for response function
00817     unsigned int num_responses = responses[i]->numResponses();
00818     Epetra_LocalMap local_response_map(num_responses, 0, comm);
00819 
00820     // Create Epetra_Vector for response function
00821     Stokhos::VectorOrthogPoly<Epetra_Vector> local_sg_g(basis,
00822               local_response_map);
00823 
00824     // Evaluate response function
00825     responses[i]->evaluateSGResponses(sg_xdot, sg_x, p, sg_p, sg_p_vals,
00826               local_sg_g);
00827 
00828     // Copy result into combined result
00829     for (int k=0; k<sg_g.size(); k++)
00830       for (unsigned int j=0; j<num_responses; j++)
00831   sg_g[k][offset+j] = local_sg_g[k][j];
00832 
00833     // Increment offset in combined result
00834     offset += num_responses;
00835   }
00836 }
00837 
00838 #endif

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