Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
00030 
00031 #include "SundanceLinearProblem.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "SundanceAssembler.hpp"
00035 #include "SundanceDiscreteFunction.hpp"
00036 #include "SundanceEquationSet.hpp"
00037 #include "SundanceZeroExpr.hpp"
00038 #include "SundanceExpr.hpp"
00039 #include "SundanceListExpr.hpp"
00040 #include "PlayaSolverState.hpp"
00041 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00042 #include "PlayaLinearOperatorImpl.hpp"
00043 #endif
00044 
00045 
00046 using namespace Sundance;
00047 using namespace Teuchos;
00048 using namespace Playa;
00049 using namespace std;
00050 
00051 
00052 static Time& lpCtorTimer() 
00053 {
00054   static RCP<Time> rtn 
00055     = TimeMonitor::getNewTimer("LinearProblem ctor"); 
00056   return *rtn;
00057 }
00058 
00059 
00060 LinearProblem::LinearProblem() 
00061   : assembler_(),
00062     A_(),
00063     rhs_()
00064 {
00065   TimeMonitor timer(lpCtorTimer());
00066 }
00067 
00068 
00069 
00070 LinearProblem::LinearProblem(const Mesh& mesh, 
00071   const Expr& eqn, 
00072   const Expr& bc,
00073   const Expr& test, 
00074   const Expr& unk, 
00075   const VectorType<double>& vecType)
00076   : assembler_(),
00077     A_(),
00078     rhs_(1),
00079     names_(1),
00080     solveDriver_(),
00081     params_()
00082 {
00083   bool partitionBCs = false;
00084   TimeMonitor timer(lpCtorTimer());
00085   Expr u = unk.flattenSpectral();
00086   Expr v = test.flattenSpectral();
00087 
00088   Array<Expr> zero(u.size());
00089   for (int i=0; i<u.size(); i++) 
00090   {
00091     Expr z = new ZeroExpr();
00092     zero[i] = z;
00093     names_[0].append(u[i].toString());
00094   }
00095 
00096   Expr u0 = new ListExpr(zero);
00097 
00098   Expr unkParams;
00099   Expr fixedParams;
00100   Array<Expr> fixedFields;
00101   Expr unkParamValues;
00102   Expr fixedParamValues;
00103   Array<Expr> fixedFieldValues;
00104 
00105 
00106   RCP<EquationSet> eqnSet 
00107     = rcp(new EquationSet(eqn, bc, tuple(v), tuple(u), tuple(u0),
00108         unkParams, unkParamValues,
00109         fixedParams, fixedParamValues,
00110         fixedFields, fixedFieldValues));
00111 
00112   assembler_ = rcp(new Assembler(mesh, eqnSet, tuple(vecType), tuple(vecType), partitionBCs));
00113 }
00114 
00115 
00116 LinearProblem::LinearProblem(const Mesh& mesh, 
00117   const Expr& eqn, 
00118   const Expr& bc,
00119   const Expr& test, 
00120   const Expr& unk, 
00121   const Expr& params, 
00122   const Expr& paramVals, 
00123   const VectorType<double>& vecType)
00124   : assembler_(),
00125     A_(),
00126     rhs_(1),
00127     names_(1),
00128     solveDriver_(),
00129     params_(params)
00130 {
00131   bool partitionBCs = false;
00132   TimeMonitor timer(lpCtorTimer());
00133   Expr u = unk.flattenSpectral();
00134   Expr v = test.flattenSpectral();
00135 
00136   Expr dumParams;
00137   Expr alpha = params.flattenSpectral();
00138   Expr alpha0 = paramVals.flattenSpectral();
00139   Array<Expr> zero(u.size());
00140   for (int i=0; i<u.size(); i++) 
00141   {
00142     Expr z = new ZeroExpr();
00143     zero[i] = z;
00144     names_[0].append(u[i].toString());
00145   }
00146 
00147   Expr u0 = new ListExpr(zero);
00148 
00149   Array<Expr> fixedFields;
00150 
00151   
00152   RCP<EquationSet> eqnSet 
00153     = rcp(new EquationSet(eqn, bc, tuple(v), tuple(u), tuple(u0),
00154         dumParams, dumParams, 
00155         alpha, alpha0,
00156         fixedFields, fixedFields));
00157 
00158   assembler_ = rcp(new Assembler(mesh, eqnSet, tuple(vecType), tuple(vecType), partitionBCs));
00159 }
00160 
00161 
00162 
00163 LinearProblem::LinearProblem(const Mesh& mesh, 
00164   const Expr& eqn, 
00165   const Expr& bc,
00166   const BlockArray& test, 
00167   const BlockArray& unk)
00168   : assembler_(),
00169     A_(),
00170     rhs_(1),
00171     names_(unk.size()),
00172     solveDriver_(),
00173     params_()
00174 {
00175   bool partitionBCs = false;
00176   TimeMonitor timer(lpCtorTimer());
00177   Array<Expr> v(test.size());  
00178   Array<Expr> u(unk.size());
00179   Array<Expr> u0(unk.size());
00180 
00181   Array<VectorType<double> > testVecType(test.size());
00182   Array<VectorType<double> > unkVecType(unk.size());
00183 
00184   for (int i=0; i<test.size(); i++)
00185   {
00186     v[i] = test[i].expr().flattenSpectral();
00187     testVecType[i] = test[i].vecType();
00188   }
00189 
00190   for (int i=0; i<unk.size(); i++)
00191   {
00192     u[i] = unk[i].expr().flattenSpectral();
00193     unkVecType[i] = unk[i].vecType();
00194     Array<Expr> zero(u[i].size());
00195     for (int j=0; j<u[i].size(); j++) 
00196     {
00197       Expr z = new ZeroExpr();
00198       zero[j] = z;
00199       names_[i].append(u[i][j].toString());
00200     }
00201     u0[i] = new ListExpr(zero);
00202 
00203   }
00204 
00205   Expr unkParams;
00206   Expr fixedParams;
00207   Array<Expr> fixedFields;
00208   Expr unkParamValues;
00209   Expr fixedParamValues;
00210   Array<Expr> fixedFieldValues;
00211   
00212   RCP<EquationSet> eqnSet 
00213     = rcp(new EquationSet(eqn, bc, v, u, u0,
00214         unkParams, unkParamValues,
00215         fixedParams, fixedParamValues,
00216         fixedFields, fixedFieldValues));
00217 
00218   assembler_ = rcp(new Assembler(mesh, eqnSet, testVecType, unkVecType, partitionBCs));
00219 }
00220 
00221 
00222 LinearProblem::LinearProblem(const Mesh& mesh, 
00223   const Expr& eqn, 
00224   const Expr& bc,
00225   const BlockArray& test, 
00226   const BlockArray& unk,
00227   const Expr& unkParams, 
00228   const Expr& unkParamVals)
00229   : assembler_(),
00230     A_(),
00231     rhs_(1),
00232     names_(unk.size()),
00233     solveDriver_(),
00234     params_(unkParams)
00235 {
00236   bool partitionBCs = false;
00237   TimeMonitor timer(lpCtorTimer());
00238   Array<Expr> v(test.size());  
00239   Array<Expr> u(unk.size());
00240   Array<Expr> u0(unk.size());
00241   Array<VectorType<double> > testVecType(test.size());
00242   Array<VectorType<double> > unkVecType(unk.size());
00243 
00244   for (int i=0; i<test.size(); i++)
00245   {
00246     v[i] = test[i].expr().flattenSpectral();
00247     testVecType[i] = test[i].vecType();
00248   }
00249 
00250   for (int i=0; i<unk.size(); i++)
00251   {
00252     u[i] = unk[i].expr().flattenSpectral();
00253     unkVecType[i] = unk[i].vecType();
00254     Array<Expr> zero(u[i].size());
00255     for (int j=0; j<u[i].size(); j++) 
00256     {
00257       Expr z = new ZeroExpr();
00258       zero[j] = z;
00259       names_[i].append(u[i][j].toString());
00260     }
00261     u0[i] = new ListExpr(zero);
00262   }
00263 
00264   Expr fixedParams;
00265   Array<Expr> fixedFields;
00266   Expr fixedParamValues;
00267   Array<Expr> fixedFieldValues;
00268   
00269   RCP<EquationSet> eqnSet 
00270     = rcp(new EquationSet(eqn, bc, v, u, u0,
00271         unkParams.flattenSpectral(), 
00272         unkParamVals.flattenSpectral(),
00273         fixedParams, fixedParamValues,
00274         fixedFields, fixedFieldValues));
00275 
00276   assembler_ = rcp(new Assembler(mesh, eqnSet, testVecType, unkVecType, partitionBCs));
00277 }
00278 
00279 LinearProblem::LinearProblem(const RCP<Assembler>& assembler)
00280   : assembler_(assembler),
00281     A_(),
00282     rhs_(1),
00283     names_(),
00284     params_()
00285 {  
00286   TimeMonitor timer(lpCtorTimer());
00287   const RCP<EquationSet>& eqn = assembler->eqnSet();
00288   names_.resize(eqn->numUnkBlocks());
00289   for (int i=0; i<eqn->numUnkBlocks(); i++)
00290   {
00291     for (int j=0; j<eqn->numUnks(i); j++) 
00292     {
00293       names_[i].append("u(" + Teuchos::toString(i) + ", "
00294         + Teuchos::toString(j) + ")");
00295     }
00296   }
00297 }
00298 
00299 /* Return the map from cells and functions to row indices */
00300 const RCP<DOFMapBase>& LinearProblem::rowMap(int blockRow) const 
00301 {return assembler_->rowMap()[blockRow];}
00302     
00303 /* Return the map from cells and functions to column indices */
00304 const RCP<DOFMapBase>& LinearProblem::colMap(int blockCol) const 
00305 {return assembler_->colMap()[blockCol];}
00306 
00307 /* Return the discrete space in which solutions live */
00308 const Array<RCP<DiscreteSpace> >& LinearProblem::solnSpace() const 
00309 {return assembler_->solutionSpace();}
00310     
00311 /* Return the set of row indices marked as 
00312  * essential boundary conditions */
00313 const RCP<Set<int> >& LinearProblem::bcRows(int blockRow) const 
00314 {return assembler_->bcRows()[blockRow];}
00315 
00316 /* Return the number of block rows in the problem  */
00317 int LinearProblem::numBlockRows() const {return assembler_->rowMap().size();}
00318 
00319 /* Return the number of block cols in the problem  */
00320 int LinearProblem::numBlockCols() const {return assembler_->colMap().size();}
00321 
00322 Array<Vector<double> > LinearProblem::getRHS() const 
00323 {
00324   Tabs tab;
00325   SUNDANCE_MSG1(assembler_->maxWatchFlagSetting("solve control"),
00326     tab << "LinearProblem::solve() building vector");
00327   assembler_->assemble(rhs_);
00328   return rhs_;
00329 }
00330 
00331 
00332 Playa::LinearOperator<double> LinearProblem::getOperator() const 
00333 {
00334   Tabs tab;
00335   SUNDANCE_MSG1(assembler_->maxWatchFlagSetting("solve control"),
00336     tab << "LinearProblem::solve() building matrix and vector");
00337   assembler_->assemble(A_, rhs_);
00338   return A_;
00339 }
00340 
00341 Expr LinearProblem::solve(const LinearSolver<double>& solver) const 
00342 {
00343   Tabs tab;
00344   int verb = assembler_->maxWatchFlagSetting("solve control");
00345   Array<Vector<double> > solnVec(rhs_.size());
00346   
00347   SUNDANCE_MSG1(verb, tab << "LinearProblem::solve() building system");
00348 
00349   assembler_->assemble(A_, rhs_);
00350   for (int i=0; i<rhs_.size(); i++)
00351     rhs_[i].scale(-1.0);
00352 
00353   SUNDANCE_MSG1(verb, tab << "LinearProblem::solve() solving system");
00354 
00355   Expr rtn;
00356 
00357   /* we're not checking the status of the solve, so failures should
00358    * be considered fatal */
00359   bool save = LinearSolveDriver::solveFailureIsFatal();
00360   LinearSolveDriver::solveFailureIsFatal() = true;
00361 
00362   solveDriver_.solve(solver, A_, rhs_, solnSpace(), names_, verb, rtn);
00363 
00364   SUNDANCE_MSG1(verb, tab << "LinearProblem::solve() done solving system");
00365 
00366   /* restore original failure-handling setting */
00367   LinearSolveDriver::solveFailureIsFatal()=save; 
00368 
00369   return rtn;
00370 }
00371 
00372 SolverState<double> LinearProblem
00373 ::solve(const LinearSolver<double>& solver,
00374   Expr& soln) const 
00375 {
00376   Tabs tab;
00377   int verb = assembler_->maxWatchFlagSetting("solve control");
00378 
00379   Array<Vector<double> > solnVec(rhs_.size());
00380   
00381   SUNDANCE_MSG1(verb, tab << "LinearProblem::solve() building system");
00382 
00383   assembler_->assemble(A_, rhs_);
00384   for (int i=0; i<rhs_.size(); i++)
00385   {
00386     rhs_[i].scale(-1.0);
00387   }
00388 
00389   SUNDANCE_MSG1(verb, tab << "solving LinearProblem");
00390   
00391   return solveDriver_.solve(solver, A_, rhs_, solnSpace(), names_, verb, soln);
00392 }
00393 
00394 
00395 Expr LinearProblem::formSolutionExpr(const Array<Vector<double> >& vec) const
00396 {
00397   int verb = assembler_->maxWatchFlagSetting("solve control");
00398   return solveDriver_.formSolutionExpr(vec, solnSpace(), names_, verb);
00399 }
00400 
00401 
00402 void LinearProblem::reAssembleProblem() const {
00403   assembler_->flushConfiguration();
00404 }
00405
SundanceLinearProblem.cpp
