Albany: a Trilinos-based PDE code: Albany_SolutionResponseFunction.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_SolutionResponseFunction.hpp"
00008 #include "Epetra_CrsMatrix.h"
00009 #include <algorithm>
00010 
00011 Albany::SolutionResponseFunction::
00012 SolutionResponseFunction(
00013   const Teuchos::RCP<Albany::Application>& application_,
00014   Teuchos::ParameterList& responseParams) : 
00015   application(application_)
00016 {
00017   // Build list of DOFs we want to keep
00018   // This should be replaced by DOF names eventually
00019   int numDOF = application->getProblem()->numEquations();
00020   if (responseParams.isType< Teuchos::Array<int> >("Keep DOF Indices")) {
00021     Teuchos::Array<int> dofs = 
00022       responseParams.get< Teuchos::Array<int> >("Keep DOF Indices");
00023     keepDOF = Teuchos::Array<int>(numDOF, 0);
00024     for (int i=0; i<dofs.size(); i++)
00025       keepDOF[dofs[i]] = 1;
00026   }
00027   else {
00028     keepDOF = Teuchos::Array<int>(numDOF, 1);
00029   }
00030 }
00031 
00032 void
00033 Albany::SolutionResponseFunction::
00034 setup()
00035 {
00036   // Build culled map and importer
00037   Teuchos::RCP<const Epetra_Map> x_map = application->getMap();
00038   culled_map = buildCulledMap(*x_map, keepDOF);
00039   importer = Teuchos::rcp(new Epetra_Import(*culled_map, *x_map));
00040 
00041   // Create graph for gradient operator -- diagonal matrix
00042   gradient_graph = 
00043     Teuchos::rcp(new Epetra_CrsGraph(Copy, *culled_map, 1, true));
00044   for (int i=0; i<culled_map->NumMyElements(); i++) {
00045     int row = culled_map->GID(i);
00046     gradient_graph->InsertGlobalIndices(row, 1, &row);
00047   }
00048   gradient_graph->FillComplete();
00049   gradient_graph->OptimizeStorage();
00050 }
00051 
00052 Albany::SolutionResponseFunction::
00053 ~SolutionResponseFunction()
00054 {
00055 }
00056 
00057 Teuchos::RCP<const Epetra_Map>
00058 Albany::SolutionResponseFunction::
00059 responseMap() const
00060 {
00061   return culled_map;
00062 }
00063 
00064 Teuchos::RCP<Epetra_Operator>
00065 Albany::SolutionResponseFunction::
00066 createGradientOp() const
00067 {
00068   Teuchos::RCP<Epetra_CrsMatrix> G =
00069     Teuchos::rcp(new Epetra_CrsMatrix(Copy, *gradient_graph));
00070   G->FillComplete();
00071   return G;
00072 }
00073 
00074 void
00075 Albany::SolutionResponseFunction::
00076 evaluateResponse(const double current_time,
00077      const Epetra_Vector* xdot,
00078      const Epetra_Vector* xdotdot,
00079      const Epetra_Vector& x,
00080      const Teuchos::Array<ParamVec>& p,
00081      Epetra_Vector& g)
00082 {
00083   cullSolution(x, g);
00084 }
00085 
00086 void
00087 Albany::SolutionResponseFunction::
00088 evaluateTangent(const double alpha, 
00089     const double beta,
00090     const double omega,
00091     const double current_time,
00092     bool sum_derivs,
00093     const Epetra_Vector* xdot,
00094     const Epetra_Vector* xdotdot,
00095     const Epetra_Vector& x,
00096     const Teuchos::Array<ParamVec>& p,
00097     ParamVec* deriv_p,
00098     const Epetra_MultiVector* Vxdot,
00099     const Epetra_MultiVector* Vxdotdot,
00100     const Epetra_MultiVector* Vx,
00101     const Epetra_MultiVector* Vp,
00102     Epetra_Vector* g,
00103     Epetra_MultiVector* gx,
00104     Epetra_MultiVector* gp)
00105 {
00106   if (g)
00107     cullSolution(x, *g);
00108 
00109   if (gx) {
00110     gx->PutScalar(0.0);
00111     if (Vx) {
00112       cullSolution(*Vx, *gx);
00113       gx->Scale(beta);
00114     }
00115   }
00116 
00117   if (gp)
00118     gp->PutScalar(0.0);
00119 }
00120 
00121 void
00122 Albany::SolutionResponseFunction::
00123 evaluateGradient(const double current_time,
00124      const Epetra_Vector* xdot,
00125      const Epetra_Vector* xdotdot,
00126      const Epetra_Vector& x,
00127      const Teuchos::Array<ParamVec>& p,
00128      ParamVec* deriv_p,
00129      Epetra_Vector* g,
00130      Epetra_Operator* dg_dx,
00131      Epetra_Operator* dg_dxdot,
00132      Epetra_Operator* dg_dxdotdot,
00133      Epetra_MultiVector* dg_dp)
00134 {
00135   if (g)
00136     cullSolution(x, *g);
00137 
00138   if (dg_dx) {
00139     Epetra_CrsMatrix *dg_dx_crs = dynamic_cast<Epetra_CrsMatrix*>(dg_dx);
00140     TEUCHOS_TEST_FOR_EXCEPT(dg_dx_crs == NULL);
00141     dg_dx_crs->PutScalar(1.0); // matrix only stores the diagonal
00142   }
00143 
00144   if (dg_dxdot) {
00145     Epetra_CrsMatrix *dg_dxdot_crs = dynamic_cast<Epetra_CrsMatrix*>(dg_dxdot);
00146     TEUCHOS_TEST_FOR_EXCEPT(dg_dxdot_crs == NULL);
00147     dg_dxdot_crs->PutScalar(0.0); // matrix only stores the diagonal
00148   }
00149   if (dg_dxdotdot) {
00150     Epetra_CrsMatrix *dg_dxdotdot_crs = dynamic_cast<Epetra_CrsMatrix*>(dg_dxdotdot);
00151     TEUCHOS_TEST_FOR_EXCEPT(dg_dxdotdot_crs == NULL);
00152     dg_dxdotdot_crs->PutScalar(0.0); // matrix only stores the diagonal
00153   }
00154 
00155   if (dg_dp)
00156     dg_dp->PutScalar(0.0);
00157 }
00158 
00159 #ifdef ALBANY_SG_MP
00160 void
00161 Albany::SolutionResponseFunction::
00162 init_sg(
00163   const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis,
00164   const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& quad,
00165   const Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> >& expansion,
00166   const Teuchos::RCP<const EpetraExt::MultiComm>& multiComm)
00167 {
00168 }
00169 
00170 void
00171 Albany::SolutionResponseFunction::
00172 evaluateSGResponse(const double curr_time,
00173        const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00174        const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00175        const Stokhos::EpetraVectorOrthogPoly& sg_x,
00176        const Teuchos::Array<ParamVec>& p,
00177        const Teuchos::Array<int>& sg_p_index,
00178        const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00179        Stokhos::EpetraVectorOrthogPoly& sg_g)
00180 {
00181   // Note that by doing the culling this way, instead of importing into
00182   // sg_g directly using a product importer, it doesn't really matter that
00183   // the product maps between sg_x and sg_g aren't consistent
00184   for (int i=0; i<sg_g.size(); i++)
00185     cullSolution(sg_x[i], sg_g[i]);
00186 }
00187 
00188 void
00189 Albany::SolutionResponseFunction::
00190 evaluateSGTangent(const double alpha, 
00191       const double beta, 
00192       const double omega, 
00193       const double current_time,
00194       bool sum_derivs,
00195       const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00196       const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00197       const Stokhos::EpetraVectorOrthogPoly& sg_x,
00198       const Teuchos::Array<ParamVec>& p,
00199       const Teuchos::Array<int>& sg_p_index,
00200       const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00201       ParamVec* deriv_p,
00202       const Epetra_MultiVector* Vx,
00203       const Epetra_MultiVector* Vxdot,
00204       const Epetra_MultiVector* Vxdotdot,
00205       const Epetra_MultiVector* Vp,
00206       Stokhos::EpetraVectorOrthogPoly* sg_g,
00207       Stokhos::EpetraMultiVectorOrthogPoly* sg_JV,
00208       Stokhos::EpetraMultiVectorOrthogPoly* sg_gp)
00209 {
00210   if (sg_g)
00211     for (int i=0; i<sg_g->size(); i++)
00212       cullSolution(sg_x[i], (*sg_g)[i]);
00213 
00214   if (sg_JV) {
00215     sg_JV->init(0.0);
00216     if (Vx) {
00217       cullSolution(*Vx, (*sg_JV)[0]);
00218       (*sg_JV)[0].Scale(beta);
00219     }
00220   }
00221 
00222   if (sg_gp)
00223     sg_gp->init(0.0);
00224 }
00225 
00226 void
00227 Albany::SolutionResponseFunction::
00228 evaluateSGGradient(const double current_time,
00229        const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00230        const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00231        const Stokhos::EpetraVectorOrthogPoly& sg_x,
00232        const Teuchos::Array<ParamVec>& p,
00233        const Teuchos::Array<int>& sg_p_index,
00234        const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00235        ParamVec* deriv_p,
00236        Stokhos::EpetraVectorOrthogPoly* sg_g,
00237        Stokhos::EpetraOperatorOrthogPoly* sg_dg_dx,
00238        Stokhos::EpetraOperatorOrthogPoly* sg_dg_dxdot,
00239        Stokhos::EpetraOperatorOrthogPoly* sg_dg_dxdotdot,
00240        Stokhos::EpetraMultiVectorOrthogPoly* sg_dg_dp)
00241 {
00242   if (sg_g)
00243     for (int i=0; i<sg_g->size(); i++)
00244       cullSolution(sg_x[i], (*sg_g)[i]);
00245 
00246   if (sg_dg_dx) {
00247     sg_dg_dx->init(0.0);
00248     Teuchos::RCP<Epetra_CrsMatrix> dg_dx_crs = 
00249       Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_dg_dx->getCoeffPtr(0),
00250               true);
00251     dg_dx_crs->PutScalar(1.0); // matrix only stores the diagonal
00252   }
00253 
00254   if (sg_dg_dxdot) {
00255     sg_dg_dxdot->init(0.0);
00256   }
00257   if (sg_dg_dxdotdot) {
00258     sg_dg_dxdotdot->init(0.0);
00259   }
00260 
00261   if (sg_dg_dp)
00262     sg_dg_dp->init(0.0);
00263 }
00264 
00265 void
00266 Albany::SolutionResponseFunction::
00267 evaluateMPResponse(const double curr_time,
00268        const Stokhos::ProductEpetraVector* mp_xdot,
00269        const Stokhos::ProductEpetraVector* mp_xdotdot,
00270        const Stokhos::ProductEpetraVector& mp_x,
00271        const Teuchos::Array<ParamVec>& p,
00272        const Teuchos::Array<int>& mp_p_index,
00273        const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00274        Stokhos::ProductEpetraVector& mp_g)
00275 {
00276   for (int i=0; i<mp_g.size(); i++)
00277     cullSolution(mp_x[i], mp_g[i]);
00278 }
00279 
00280 
00281 void
00282 Albany::SolutionResponseFunction::
00283 evaluateMPTangent(const double alpha, 
00284       const double beta, 
00285       const double omega, 
00286       const double current_time,
00287       bool sum_derivs,
00288       const Stokhos::ProductEpetraVector* mp_xdot,
00289       const Stokhos::ProductEpetraVector* mp_xdotdot,
00290       const Stokhos::ProductEpetraVector& mp_x,
00291       const Teuchos::Array<ParamVec>& p,
00292       const Teuchos::Array<int>& mp_p_index,
00293       const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00294       ParamVec* deriv_p,
00295       const Epetra_MultiVector* Vx,
00296       const Epetra_MultiVector* Vxdot,
00297       const Epetra_MultiVector* Vxdotdot,
00298       const Epetra_MultiVector* Vp,
00299       Stokhos::ProductEpetraVector* mp_g,
00300       Stokhos::ProductEpetraMultiVector* mp_JV,
00301       Stokhos::ProductEpetraMultiVector* mp_gp)
00302 {
00303   if (mp_g)
00304     for (int i=0; i<mp_g->size(); i++)
00305       cullSolution(mp_x[i], (*mp_g)[i]);
00306 
00307   if (mp_JV) {
00308     mp_JV->init(0.0);
00309     if (Vx) {
00310       for (int i=0; i<mp_JV->size(); i++) {
00311   cullSolution(*Vx, (*mp_JV)[i]);
00312   (*mp_JV)[i].Scale(beta);
00313       }
00314     }
00315   }
00316 
00317   if (mp_gp)
00318     mp_gp->init(0.0);
00319 }
00320 
00321 void
00322 Albany::SolutionResponseFunction::
00323 evaluateMPGradient(const double current_time,
00324        const Stokhos::ProductEpetraVector* mp_xdot,
00325        const Stokhos::ProductEpetraVector* mp_xdotdot,
00326        const Stokhos::ProductEpetraVector& mp_x,
00327        const Teuchos::Array<ParamVec>& p,
00328        const Teuchos::Array<int>& mp_p_index,
00329        const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00330        ParamVec* deriv_p,
00331        Stokhos::ProductEpetraVector* mp_g,
00332        Stokhos::ProductEpetraOperator* mp_dg_dx,
00333        Stokhos::ProductEpetraOperator* mp_dg_dxdot,
00334        Stokhos::ProductEpetraOperator* mp_dg_dxdotdot,
00335        Stokhos::ProductEpetraMultiVector* mp_dg_dp)
00336 {
00337   if (mp_g)
00338     for (int i=0; i<mp_g->size(); i++)
00339       cullSolution(mp_x[i], (*mp_g)[i]);
00340 
00341   if (mp_dg_dx) {
00342     for (int i=0; i<mp_dg_dx->size(); i++) {
00343       Teuchos::RCP<Epetra_CrsMatrix> dg_dx_crs = 
00344   Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(mp_dg_dx->getCoeffPtr(i),
00345                 true);
00346       dg_dx_crs->PutScalar(1.0); // matrix only stores the diagonal
00347     }
00348   }
00349 
00350   if (mp_dg_dxdot) {
00351     mp_dg_dxdot->init(0.0);
00352   }
00353   if (mp_dg_dxdotdot) {
00354     mp_dg_dxdotdot->init(0.0);
00355   }
00356 
00357   if (mp_dg_dp)
00358     mp_dg_dp->init(0.0);
00359 }
00360 #endif //ALBANY_SG_MP
00361 
00362 Teuchos::RCP<Epetra_Map> 
00363 Albany::SolutionResponseFunction::
00364 buildCulledMap(const Epetra_Map& x_map, 
00365          const Teuchos::Array<int>& keepDOF) const 
00366 {
00367   int numKeepDOF = std::accumulate(keepDOF.begin(), keepDOF.end(), 0);
00368   int Neqns = keepDOF.size();
00369   int N = x_map.NumMyElements(); // x_map is map for solution vector
00370 
00371   TEUCHOS_ASSERT( !(N % Neqns) ); // Assume that all the equations for
00372                                   // a given node are on the assigned
00373                                   // processor. I.e. need to ensure
00374                                   // that N is exactly Neqns-divisible
00375 
00376   int nnodes = N / Neqns;          // number of fem nodes
00377   int N_new = nnodes * numKeepDOF; // length of local x_new   
00378 
00379   int *gids = x_map.MyGlobalElements(); // Fill local x_map into gids array
00380   Teuchos::Array<int> gids_new(N_new);
00381   int idx = 0;
00382   for ( int inode = 0; inode < N/Neqns ; ++inode) // For every node 
00383     for ( int ieqn = 0; ieqn < Neqns; ++ieqn )  // Check every dof on the node
00384       if ( keepDOF[ieqn] == 1 )  // then want to keep this dof
00385   gids_new[idx++] = gids[(inode*Neqns)+ieqn];
00386   // end cull
00387   
00388   Teuchos::RCP<Epetra_Map> x_new_map = 
00389     Teuchos::rcp( new Epetra_Map( -1, N_new, &gids_new[0], 0, x_map.Comm() ) );
00390   
00391   return x_new_map;
00392 }
00393 
00394 void
00395 Albany::SolutionResponseFunction::
00396 cullSolution(const Epetra_MultiVector& x, Epetra_MultiVector& x_culled) const
00397 {
00398   x_culled.Import(x, *importer, Insert);
00399 }