Albany: a Trilinos-based PDE code: Albany_KLResponseFunction.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_KLResponseFunction.hpp"
00008 #include "Stokhos_PCEAnasaziKL.hpp"
00009 #include "Teuchos_Array.hpp"
00010 #include "Teuchos_VerboseObject.hpp"
00011 
00012 Albany::KLResponseFunction::
00013 KLResponseFunction(
00014   const Teuchos::RCP<Albany::AbstractResponseFunction>& response_,
00015   Teuchos::ParameterList& responseParams) :
00016   response(response_),
00017   responseParams(responseParams),
00018   out(Teuchos::VerboseObjectBase::getDefaultOStream())
00019 {
00020   num_kl = responseParams.get("Number of KL Terms", 5);
00021 }
00022 
00023 Albany::KLResponseFunction::
00024 ~KLResponseFunction()
00025 {
00026 }
00027 
00028 Teuchos::RCP<const Epetra_Map>
00029 Albany::KLResponseFunction::
00030 responseMap() const
00031 {
00032   return response->responseMap();
00033 }
00034 
00035 Teuchos::RCP<Epetra_Operator> 
00036 Albany::KLResponseFunction::
00037 createGradientOp() const
00038 {
00039   return response->createGradientOp();
00040 }
00041 
00042 bool
00043 Albany::KLResponseFunction::
00044 isScalarResponse() const
00045 {
00046   return response->isScalarResponse();
00047 }
00048 
00049 void
00050 Albany::KLResponseFunction::
00051 evaluateResponse(const double current_time,
00052      const Epetra_Vector* xdot,
00053      const Epetra_Vector* xdotdot,
00054      const Epetra_Vector& x,
00055      const Teuchos::Array<ParamVec>& p,
00056      Epetra_Vector& g)
00057 {
00058   response->evaluateResponse(current_time, xdot, xdotdot, x, p, g);
00059 }
00060 
00061 void
00062 Albany::KLResponseFunction::
00063 evaluateTangent(const double alpha, 
00064     const double beta,
00065     const double omega,
00066     const double current_time,
00067     bool sum_derivs,
00068     const Epetra_Vector* xdot,
00069     const Epetra_Vector* xdotdot,
00070     const Epetra_Vector& x,
00071     const Teuchos::Array<ParamVec>& p,
00072     ParamVec* deriv_p,
00073     const Epetra_MultiVector* Vxdot,
00074     const Epetra_MultiVector* Vxdotdot,
00075     const Epetra_MultiVector* Vx,
00076     const Epetra_MultiVector* Vp,
00077     Epetra_Vector* g,
00078     Epetra_MultiVector* gx,
00079     Epetra_MultiVector* gp)
00080 {
00081   response->evaluateTangent(alpha, beta, omega, current_time, sum_derivs, 
00082           xdot, xdotdot, x, p, deriv_p, Vxdot, Vxdotdot, Vx, Vp,
00083           g, gx, gp);
00084 }
00085 
00086 void
00087 Albany::KLResponseFunction::
00088 evaluateDerivative(const double current_time,
00089        const Epetra_Vector* xdot,
00090        const Epetra_Vector* xdotdot,
00091        const Epetra_Vector& x,
00092        const Teuchos::Array<ParamVec>& p,
00093        ParamVec* deriv_p,
00094        Epetra_Vector* g,
00095        const EpetraExt::ModelEvaluator::Derivative& dg_dx,
00096        const EpetraExt::ModelEvaluator::Derivative& dg_dxdot,
00097        const EpetraExt::ModelEvaluator::Derivative& dg_dxdotdot,
00098        const EpetraExt::ModelEvaluator::Derivative& dg_dp)
00099 {
00100   response->evaluateDerivative(current_time, xdot, xdotdot, x, p, deriv_p, 
00101              g, dg_dx, dg_dxdot, dg_dxdotdot, dg_dp);
00102 }
00103 
00104 #ifdef ALBANY_SG_MP
00105 void
00106 Albany::KLResponseFunction::
00107 init_sg(
00108   const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis,
00109   const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& quad,
00110   const Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> >& expansion,
00111   const Teuchos::RCP<const EpetraExt::MultiComm>& multiComm)
00112 {
00113 }
00114 
00115 void
00116 Albany::KLResponseFunction::
00117 evaluateSGResponse(const double curr_time,
00118        const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00119        const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00120        const Stokhos::EpetraVectorOrthogPoly& sg_x,
00121        const Teuchos::Array<ParamVec>& p,
00122        const Teuchos::Array<int>& sg_p_index,
00123        const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00124        Stokhos::EpetraVectorOrthogPoly& sg_g)
00125 {
00126   // Compute base response
00127   response->evaluateSGResponse(curr_time, sg_xdot, sg_xdotdot, sg_x, p, sg_p_index,
00128              sg_p_vals, sg_g);
00129 
00130   // Compute KL of that response
00131   Teuchos::Array<double> evals;
00132   Teuchos::RCP<Epetra_MultiVector> evecs;
00133   bool success = computeKL(sg_g, num_kl, evals, evecs);
00134 
00135   if (!success) 
00136     *out << "KL Eigensolver did not converge!" << std::endl;
00137   
00138   *out << "Eigenvalues = " << std::endl;
00139   for (int i=0; i<num_kl; i++)
00140     *out << evals[i] << std::endl;
00141 }
00142 
00143 void
00144 Albany::KLResponseFunction::
00145 evaluateSGTangent(const double alpha, 
00146       const double beta, 
00147       const double omega, 
00148       const double current_time,
00149       bool sum_derivs,
00150       const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00151       const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00152       const Stokhos::EpetraVectorOrthogPoly& sg_x,
00153       const Teuchos::Array<ParamVec>& p,
00154       const Teuchos::Array<int>& sg_p_index,
00155       const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00156       ParamVec* deriv_p,
00157       const Epetra_MultiVector* Vx,
00158       const Epetra_MultiVector* Vxdot,
00159       const Epetra_MultiVector* Vxdotdot,
00160       const Epetra_MultiVector* Vp,
00161       Stokhos::EpetraVectorOrthogPoly* sg_g,
00162       Stokhos::EpetraMultiVectorOrthogPoly* sg_JV,
00163       Stokhos::EpetraMultiVectorOrthogPoly* sg_gp)
00164 {
00165   // Currently we only know how to do KL on response
00166   response->evaluateSGTangent(alpha, beta, omega, current_time, sum_derivs, 
00167             sg_xdot, sg_xdotdot, sg_x, p, sg_p_index, sg_p_vals, deriv_p, 
00168             Vxdot, Vxdotdot, Vx, Vp,
00169             sg_g, sg_JV, sg_gp);
00170   if (sg_g)
00171     evaluateSGResponse(current_time, sg_xdot, sg_xdotdot, sg_x, p, sg_p_index, sg_p_vals, 
00172            *sg_g);
00173 }
00174 
00175 void
00176 Albany::KLResponseFunction::
00177 evaluateSGDerivative(
00178   const double current_time,
00179   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00180   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00181   const Stokhos::EpetraVectorOrthogPoly& sg_x,
00182   const Teuchos::Array<ParamVec>& p,
00183   const Teuchos::Array<int>& sg_p_index,
00184   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00185   ParamVec* deriv_p,
00186   Stokhos::EpetraVectorOrthogPoly* sg_g,
00187   const EpetraExt::ModelEvaluator::SGDerivative& sg_dg_dx,
00188   const EpetraExt::ModelEvaluator::SGDerivative& sg_dg_dxdot,
00189   const EpetraExt::ModelEvaluator::SGDerivative& sg_dg_dxdotdot,
00190   const EpetraExt::ModelEvaluator::SGDerivative& sg_dg_dp)
00191 {
00192   // Currently we only know how to do KL on response
00193   response->evaluateSGDerivative(current_time, sg_xdot, sg_xdotdot, sg_x, p, sg_p_index, 
00194          sg_p_vals, deriv_p, 
00195          sg_g, sg_dg_dx, sg_dg_dxdot, sg_dg_dxdotdot, sg_dg_dp);
00196 
00197   if (sg_g)
00198     evaluateSGResponse(current_time, sg_xdot, sg_xdotdot, sg_x, p, sg_p_index, sg_p_vals, 
00199            *sg_g);
00200 }
00201 
00202 void
00203 Albany::KLResponseFunction::
00204 evaluateMPResponse(const double curr_time,
00205        const Stokhos::ProductEpetraVector* mp_xdot,
00206        const Stokhos::ProductEpetraVector* mp_xdotdot,
00207        const Stokhos::ProductEpetraVector& mp_x,
00208        const Teuchos::Array<ParamVec>& p,
00209        const Teuchos::Array<int>& mp_p_index,
00210        const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00211        Stokhos::ProductEpetraVector& mp_g)
00212 {
00213   response->evaluateMPResponse(curr_time, mp_xdot, mp_xdotdot, mp_x, p, mp_p_index, 
00214              mp_p_vals, mp_g);
00215 }
00216 
00217 
00218 void
00219 Albany::KLResponseFunction::
00220 evaluateMPTangent(const double alpha, 
00221       const double beta, 
00222       const double omega, 
00223       const double current_time,
00224       bool sum_derivs,
00225       const Stokhos::ProductEpetraVector* mp_xdot,
00226       const Stokhos::ProductEpetraVector* mp_xdotdot,
00227       const Stokhos::ProductEpetraVector& mp_x,
00228       const Teuchos::Array<ParamVec>& p,
00229       const Teuchos::Array<int>& mp_p_index,
00230       const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00231       ParamVec* deriv_p,
00232       const Epetra_MultiVector* Vx,
00233       const Epetra_MultiVector* Vxdot,
00234       const Epetra_MultiVector* Vxdotdot,
00235       const Epetra_MultiVector* Vp,
00236       Stokhos::ProductEpetraVector* mp_g,
00237       Stokhos::ProductEpetraMultiVector* mp_JV,
00238       Stokhos::ProductEpetraMultiVector* mp_gp)
00239 {
00240   response->evaluateMPTangent(alpha, beta, omega, current_time, sum_derivs, 
00241             mp_xdot, mp_xdotdot, mp_x, p, mp_p_index, mp_p_vals, deriv_p, 
00242             Vxdot, Vxdotdot, Vx, Vp,
00243             mp_g, mp_JV, mp_gp);
00244 }
00245 
00246 void
00247 Albany::KLResponseFunction::
00248 evaluateMPDerivative(
00249   const double current_time,
00250   const Stokhos::ProductEpetraVector* mp_xdot,
00251   const Stokhos::ProductEpetraVector* mp_xdotdot,
00252   const Stokhos::ProductEpetraVector& mp_x,
00253   const Teuchos::Array<ParamVec>& p,
00254   const Teuchos::Array<int>& mp_p_index,
00255   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00256   ParamVec* deriv_p,
00257   Stokhos::ProductEpetraVector* mp_g,
00258   const EpetraExt::ModelEvaluator::MPDerivative& mp_dg_dx,
00259   const EpetraExt::ModelEvaluator::MPDerivative& mp_dg_dxdot,
00260   const EpetraExt::ModelEvaluator::MPDerivative& mp_dg_dxdotdot,
00261   const EpetraExt::ModelEvaluator::MPDerivative& mp_dg_dp)
00262 {
00263   response->evaluateMPDerivative(current_time, mp_xdot, mp_xdotdot, mp_x, p, mp_p_index, 
00264          mp_p_vals, deriv_p, 
00265          mp_g, mp_dg_dx, mp_dg_dxdot, mp_dg_dxdotdot, mp_dg_dp);
00266 }
00267 #endif //ALBANY_SG_MP
00268 
00269 bool
00270 Albany::KLResponseFunction::
00271 computeKL(const Stokhos::EpetraVectorOrthogPoly& sg_u,
00272     const int NumKL,
00273     Teuchos::Array<double>& evals,
00274     Teuchos::RCP<Epetra_MultiVector>& evecs)
00275 {
00276   Teuchos::RCP<const EpetraExt::BlockVector> X_ov = sg_u.getBlockVector();
00277     //App_sg->get_sg_model()->import_solution( *(sg_u->getBlockVector()) );
00278   //Teuchos::RCP<const EpetraExt::BlockVector> cX_ov = X_ov;
00279 
00280   // pceKL is object with member functions that explicitly call anasazi
00281   Stokhos::PCEAnasaziKL pceKL(X_ov, *(sg_u.basis()), NumKL);
00282 
00283   // Set parameters for anasazi
00284   Teuchos::ParameterList& anasazi_params = 
00285     responseParams.sublist("Anasazi");
00286   Teuchos::ParameterList default_params = pceKL.getDefaultParams();
00287   anasazi_params.setParametersNotAlreadySet(default_params);
00288 
00289   // Self explanatory
00290   bool result = pceKL.computeKL(anasazi_params);
00291    
00292   // Retrieve evals/evectors into return argument slots...
00293   evals = pceKL.getEigenvalues();
00294   evecs = pceKL.getEigenvectors();
00295 
00296   return result;
00297 }