Albany: a Trilinos-based PDE code: LocalNonlinearSolver

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 namespace LCM
00008 {
00009 
00010 template<typename EvalT, typename Traits>
00011 LocalNonlinearSolver_Base<EvalT, Traits>::LocalNonlinearSolver_Base() :
00012     lapack()
00013 {
00014 }
00015 
00016 // -----------------------------------------------------------------------------
00017 // Specializations
00018 // -----------------------------------------------------------------------------
00019 
00020 // -----------------------------------------------------------------------------
00021 // Residual
00022 // -----------------------------------------------------------------------------
00023 template<typename Traits>
00024 LocalNonlinearSolver<PHAL::AlbanyTraits::Residual, Traits>::LocalNonlinearSolver() :
00025     LocalNonlinearSolver_Base<PHAL::AlbanyTraits::Residual, Traits>()
00026 {
00027 }
00028 
00029 template<typename Traits>
00030 void
00031 inline
00032 LocalNonlinearSolver<PHAL::AlbanyTraits::Residual, Traits>::
00033 solve(
00034     std::vector<ScalarT> & A,
00035     std::vector<ScalarT> & X,
00036     std::vector<ScalarT> & B)
00037 {
00038   // system size
00039   int numLocalVars = B.size();
00040 
00041   // data for the LAPACK call below
00042   int info(0);
00043   std::vector<int> IPIV(numLocalVars);
00044 
00045   // call LAPACK
00046   this->lapack.GESV(numLocalVars, 1, &A[0], numLocalVars, &IPIV[0], &B[0],
00047       numLocalVars, &info);
00048 
00049   // increment the solution
00050   for (int i(0); i < numLocalVars; ++i)
00051     X[i] -= B[i];
00052 }
00053 
00054 template<typename Traits>
00055 void
00056 LocalNonlinearSolver<PHAL::AlbanyTraits::Residual, Traits>::
00057 computeFadInfo(
00058     std::vector<ScalarT> & A,
00059     std::vector<ScalarT> & X,
00060     std::vector<ScalarT> & B)
00061 {
00062   // no-op
00063 }
00064 
00065 // -----------------------------------------------------------------------------
00066 // Jacobian
00067 // -----------------------------------------------------------------------------
00068 template<typename Traits>
00069 LocalNonlinearSolver<PHAL::AlbanyTraits::Jacobian, Traits>::LocalNonlinearSolver() :
00070     LocalNonlinearSolver_Base<PHAL::AlbanyTraits::Jacobian, Traits>()
00071 {
00072 }
00073 
00074 template<typename Traits>
00075 void
00076 LocalNonlinearSolver<PHAL::AlbanyTraits::Jacobian, Traits>::
00077 solve(
00078     std::vector<ScalarT> & A,
00079     std::vector<ScalarT> & X,
00080     std::vector<ScalarT> & B)
00081 {
00082   // system size
00083   int numLocalVars = B.size();
00084 
00085   // data for the LAPACK call below
00086   int info(0);
00087   std::vector<int> IPIV(numLocalVars);
00088 
00089   // fill B and dBdX
00090   std::vector<RealType> F(numLocalVars);
00091   std::vector<RealType> dFdX(numLocalVars * numLocalVars);
00092   for (int i(0); i < numLocalVars; ++i)
00093       {
00094     F[i] = B[i].val();
00095     for (int j(0); j < numLocalVars; ++j)
00096         {
00097       dFdX[i + numLocalVars * j] = A[i + numLocalVars * j].val();
00098     }
00099   }
00100 
00101   // call LAPACK
00102   this->lapack.GESV(numLocalVars, 1, &dFdX[0], numLocalVars, &IPIV[0], &F[0],
00103       numLocalVars, &info);
00104 
00105   // increment the solution
00106   for (int i(0); i < numLocalVars; ++i)
00107     X[i].val() -= F[i];
00108 
00109 }
00110 
00111 template<typename Traits>
00112 void
00113 LocalNonlinearSolver<PHAL::AlbanyTraits::Jacobian, Traits>::
00114 computeFadInfo(
00115     std::vector<ScalarT> & A,
00116     std::vector<ScalarT> & X,
00117     std::vector<ScalarT> & B)
00118 {
00119   // local system size
00120   int numLocalVars = B.size();
00121   int numGlobalVars = B[0].size();
00122   TEUCHOS_TEST_FOR_EXCEPTION(numGlobalVars == 0, std::logic_error,
00123       "In LocalNonlinearSolver<Jacobian> the numGLobalVars is zero where it should be positive\n");
00124 
00125   // data for the LAPACK call below
00126   int info(0);
00127   std::vector<int> IPIV(numLocalVars);
00128 
00129   // extract sensitivities of objective function(s) wrt p
00130   std::vector<RealType> dBdP(numLocalVars * numGlobalVars);
00131   for (int i(0); i < numLocalVars; ++i) {
00132     for (int j(0); j < numGlobalVars; ++j) {
00133       dBdP[i + numLocalVars * j] = B[i].dx(j);
00134     }
00135   }
00136 
00137   // extract the jacobian
00138   std::vector<RealType> dBdX(A.size());
00139   for (int i(0); i < numLocalVars; ++i) {
00140     for (int j(0); j < numLocalVars; ++j) {
00141       dBdX[i + numLocalVars * j] = A[i + numLocalVars * j].val();
00142     }
00143   }
00144   // call LAPACK to simultaneously solve for all dXdP
00145   this->lapack.GESV(numLocalVars, numGlobalVars, &dBdX[0], numLocalVars,
00146       &IPIV[0], &dBdP[0], numLocalVars, &info);
00147 
00148   // unpack into globalX (recall that LAPACK stores dXdP in dBdP)
00149   for (int i(0); i < numLocalVars; ++i)
00150       {
00151     X[i].resize(numGlobalVars);
00152     for (int j(0); j < numGlobalVars; ++j)
00153         {
00154       X[i].fastAccessDx(j) = -dBdP[i + numLocalVars * j];
00155     }
00156   }
00157 }
00158 
00159 // -----------------------------------------------------------------------------
00160 // Tangent
00161 // -----------------------------------------------------------------------------
00162 template<typename Traits>
00163 LocalNonlinearSolver<PHAL::AlbanyTraits::Tangent, Traits>::LocalNonlinearSolver() :
00164     LocalNonlinearSolver_Base<PHAL::AlbanyTraits::Tangent, Traits>()
00165 {
00166 }
00167 
00168 template<typename Traits>
00169 void
00170 LocalNonlinearSolver<PHAL::AlbanyTraits::Tangent, Traits>::
00171 solve(
00172     std::vector<ScalarT> & A,
00173     std::vector<ScalarT> & X,
00174     std::vector<ScalarT> & B)
00175 {
00176   // system size
00177   int numLocalVars = B.size();
00178 
00179   // data for the LAPACK call below
00180   int info(0);
00181   std::vector<int> IPIV(numLocalVars);
00182 
00183   // fill B and dBdX
00184   std::vector<RealType> F(numLocalVars);
00185   std::vector<RealType> dFdX(numLocalVars * numLocalVars);
00186   for (int i(0); i < numLocalVars; ++i)
00187       {
00188     F[i] = B[i].val();
00189     for (int j(0); j < numLocalVars; ++j)
00190         {
00191       dFdX[i + numLocalVars * j] = A[i + numLocalVars * j].val();
00192     }
00193   }
00194 
00195   // call LAPACK
00196   this->lapack.GESV(numLocalVars, 1, &dFdX[0], numLocalVars, &IPIV[0], &F[0],
00197       numLocalVars, &info);
00198 
00199   // increment the solution
00200   for (int i(0); i < numLocalVars; ++i)
00201     X[i].val() -= F[i];
00202 }
00203 
00204 template<typename Traits>
00205 void
00206 LocalNonlinearSolver<PHAL::AlbanyTraits::Tangent, Traits>::
00207 computeFadInfo(
00208     std::vector<ScalarT> & A,
00209     std::vector<ScalarT> & X,
00210     std::vector<ScalarT> & B)
00211 {
00212   // local system size
00213   int numLocalVars = B.size();
00214   int numGlobalVars = B[0].size();
00215   TEUCHOS_TEST_FOR_EXCEPTION(numGlobalVars == 0, std::logic_error,
00216       "In LocalNonlinearSolver<Tangent, Traits> the numGLobalVars is zero where it should be positive\n");
00217 
00218   // data for the LAPACK call below
00219   int info(0);
00220   std::vector<int> IPIV(numLocalVars);
00221 
00222   // extract sensitivites of objective function(s) wrt p
00223   std::vector<RealType> dBdP(numLocalVars * numGlobalVars);
00224   for (int i(0); i < numLocalVars; ++i) {
00225     for (int j(0); j < numGlobalVars; ++j) {
00226       dBdP[i + numLocalVars * j] = B[i].dx(j);
00227     }
00228   }
00229 
00230   // extract the jacobian
00231   std::vector<RealType> dBdX(A.size());
00232   for (int i(0); i < numLocalVars; ++i) {
00233     for (int j(0); j < numLocalVars; ++j) {
00234       dBdX[i + numLocalVars * j] = A[i + numLocalVars * j].val();
00235     }
00236   }
00237 
00238   // call LAPACK to simultaneously solve for all dXdP
00239   this->lapack.GESV(numLocalVars, numGlobalVars, &dBdX[0], numLocalVars,
00240       &IPIV[0], &dBdP[0], numLocalVars, &info);
00241 
00242   // unpack into globalX (recall that LAPACK stores dXdP in dBdP)
00243   for (int i(0); i < numLocalVars; ++i)
00244       {
00245     X[i].resize(numGlobalVars);
00246     for (int j(0); j < numGlobalVars; ++j)
00247         {
00248       X[i].fastAccessDx(j) = -dBdP[i + numLocalVars * j];
00249     }
00250   }
00251 }
00252 
00253 // -----------------------------------------------------------------------------
00254 // Stochastic Galerkin Residual
00255 // -----------------------------------------------------------------------------
00256 #ifdef ALBANY_SG_MP
00257 template<typename Traits>
00258 LocalNonlinearSolver<PHAL::AlbanyTraits::SGResidual, Traits>::LocalNonlinearSolver():
00259 LocalNonlinearSolver_Base<PHAL::AlbanyTraits::SGResidual, Traits>()
00260 {}
00261 
00262 template<typename Traits>
00263 void
00264 LocalNonlinearSolver<PHAL::AlbanyTraits::SGResidual, Traits>::
00265 solve(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00266 {
00267   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Stochastic Galerkin types yet\n");
00268 }
00269 
00270 template<typename Traits>
00271 void
00272 LocalNonlinearSolver<PHAL::AlbanyTraits::SGResidual, Traits>::
00273 computeFadInfo(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00274 {
00275   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Stochastic Galerkin types yet\n");
00276 }
00277 
00278 // ---------------------------------------------------------------------
00279 // Stochastic Galerkin Jacobian
00280 // ---------------------------------------------------------------------
00281 template<typename Traits>
00282 LocalNonlinearSolver<PHAL::AlbanyTraits::SGJacobian, Traits>::LocalNonlinearSolver():
00283 LocalNonlinearSolver_Base<PHAL::AlbanyTraits::SGJacobian, Traits>()
00284 {}
00285 
00286 template<typename Traits>
00287 void
00288 LocalNonlinearSolver<PHAL::AlbanyTraits::SGJacobian, Traits>::
00289 solve(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00290 {
00291   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Stochastic Galerkin types yet\n");
00292 }
00293 
00294 template<typename Traits>
00295 void
00296 LocalNonlinearSolver<PHAL::AlbanyTraits::SGJacobian, Traits>::
00297 computeFadInfo(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00298 {
00299   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Stochastic Galerkin types yet\n");
00300 }
00301 
00302 // -----------------------------------------------------------------------------
00303 // Stochastic Galerkin Tangent
00304 // -----------------------------------------------------------------------------
00305 template<typename Traits>
00306 LocalNonlinearSolver<PHAL::AlbanyTraits::SGTangent, Traits>::LocalNonlinearSolver():
00307 LocalNonlinearSolver_Base<PHAL::AlbanyTraits::SGTangent, Traits>()
00308 {}
00309 
00310 template<typename Traits>
00311 void
00312 LocalNonlinearSolver<PHAL::AlbanyTraits::SGTangent, Traits>::
00313 solve(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00314 {
00315   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Stochastic Galerkin types yet\n");
00316 }
00317 
00318 template<typename Traits>
00319 void
00320 LocalNonlinearSolver<PHAL::AlbanyTraits::SGTangent, Traits>::
00321 computeFadInfo(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00322 {
00323   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Stochastic Galerkin types yet\n");
00324 }
00325 
00326 // -----------------------------------------------------------------------------
00327 // Multi-Point Residual
00328 // -----------------------------------------------------------------------------
00329 template<typename Traits>
00330 LocalNonlinearSolver<PHAL::AlbanyTraits::MPResidual, Traits>::LocalNonlinearSolver():
00331 LocalNonlinearSolver_Base<PHAL::AlbanyTraits::MPResidual, Traits>()
00332 {}
00333 
00334 template<typename Traits>
00335 void
00336 LocalNonlinearSolver<PHAL::AlbanyTraits::MPResidual, Traits>::
00337 solve(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00338 {
00339   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Multi-Point types yet\n");
00340 }
00341 
00342 template<typename Traits>
00343 void
00344 LocalNonlinearSolver<PHAL::AlbanyTraits::MPResidual, Traits>::
00345 computeFadInfo(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00346 {
00347   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Multi-Point types yet\n");
00348 }
00349 
00350 // -----------------------------------------------------------------------------
00351 // Multi-Point Jacobian
00352 // -----------------------------------------------------------------------------
00353 template<typename Traits>
00354 LocalNonlinearSolver<PHAL::AlbanyTraits::MPJacobian, Traits>::LocalNonlinearSolver():
00355 LocalNonlinearSolver_Base<PHAL::AlbanyTraits::MPJacobian, Traits>()
00356 {}
00357 
00358 template<typename Traits>
00359 void
00360 LocalNonlinearSolver<PHAL::AlbanyTraits::MPJacobian, Traits>::
00361 solve(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00362 {
00363   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Multi-Point types yet\n");
00364 }
00365 
00366 template<typename Traits>
00367 void
00368 LocalNonlinearSolver<PHAL::AlbanyTraits::MPJacobian, Traits>::
00369 computeFadInfo(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00370 {
00371   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Multi-Point types yet\n");
00372 }
00373 
00374 // -----------------------------------------------------------------------------
00375 // Multi-Point Tangent
00376 // -----------------------------------------------------------------------------
00377 template<typename Traits>
00378 LocalNonlinearSolver<PHAL::AlbanyTraits::MPTangent, Traits>::LocalNonlinearSolver():
00379 LocalNonlinearSolver_Base<PHAL::AlbanyTraits::MPTangent, Traits>()
00380 {}
00381 
00382 template<typename Traits>
00383 void
00384 LocalNonlinearSolver<PHAL::AlbanyTraits::MPTangent, Traits>::
00385 solve(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00386 {
00387   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Multi-Point types yet\n");
00388 }
00389 
00390 template<typename Traits>
00391 void
00392 LocalNonlinearSolver<PHAL::AlbanyTraits::MPTangent, Traits>::
00393 computeFadInfo(std::vector<ScalarT> & A, std::vector<ScalarT> & X, std::vector<ScalarT> & B)
00394 {
00395   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"LocalNonlinearSolver has not been implemented for Multi-Point types yet\n");
00396 }
00397 #endif //ALBANY_SG_MP
00398 // -----------------------------------------------------------------------------
00399 }
00400
LocalNonlinearSolver_Def.hpp
