00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
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
00300 const RCP<DOFMapBase>& LinearProblem::rowMap(int blockRow) const
00301 {return assembler_->rowMap()[blockRow];}
00302
00303
00304 const RCP<DOFMapBase>& LinearProblem::colMap(int blockCol) const
00305 {return assembler_->colMap()[blockCol];}
00306
00307
00308 const Array<RCP<DiscreteSpace> >& LinearProblem::solnSpace() const
00309 {return assembler_->solutionSpace();}
00310
00311
00312
00313 const RCP<Set<int> >& LinearProblem::bcRows(int blockRow) const
00314 {return assembler_->bcRows()[blockRow];}
00315
00316
00317 int LinearProblem::numBlockRows() const {return assembler_->rowMap().size();}
00318
00319
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
00358
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
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