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 "SundanceNLOp.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "SundanceAssembler.hpp"
00035 #include "SundanceDiscreteFunction.hpp"
00036 #include "SundanceEquationSet.hpp"
00037 #include "SundanceLinearSolveDriver.hpp"
00038 #include "PlayaLinearSolverDecl.hpp"
00039
00040 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00041 #include "PlayaLinearOperatorImpl.hpp"
00042 #include "PlayaLinearSolverImpl.hpp"
00043 #include "PlayaLinearCombinationImpl.hpp"
00044 #endif
00045
00046 using namespace Sundance;
00047 using namespace Sundance;
00048 using namespace Sundance;
00049 using namespace Sundance;
00050 using namespace Sundance;
00051 using namespace Sundance;
00052 using namespace Sundance;
00053 using namespace Teuchos;
00054 using namespace std;
00055 using namespace Playa;
00056
00057
00058 static Time& nlpCtorTimer()
00059 {
00060 static RCP<Time> rtn
00061 = TimeMonitor::getNewTimer("NLOp ctor");
00062 return *rtn;
00063 }
00064
00065
00066 NLOp::NLOp()
00067 : NonlinearOperatorBase<double>(),
00068 assembler_(),
00069 u0_(),
00070 params_()
00071
00072 {
00073 TimeMonitor timer(nlpCtorTimer());
00074 }
00075
00076
00077 NLOp::NLOp(const Mesh& mesh,
00078 const Expr& eqn,
00079 const Expr& bc,
00080 const Expr& test,
00081 const Expr& unk,
00082 const Expr& u0,
00083 const VectorType<double>& vecType,
00084 bool partitionBCs)
00085 : NonlinearOperatorBase<double>(),
00086 assembler_(),
00087 u0_(u0),
00088 params_()
00089 {
00090 TimeMonitor timer(nlpCtorTimer());
00091
00092 Expr unkParams;
00093 Expr fixedParams;
00094 Expr fixedFields;
00095 Expr unkParamValues;
00096 Expr fixedParamValues;
00097 Expr fixedFieldValues;
00098
00099 RCP<EquationSet> eqnSet
00100 = rcp(new EquationSet(eqn, bc, tuple(test.flattenSpectral()), tuple(unk.flattenSpectral()), tuple(u0),
00101 unkParams, unkParamValues,
00102 fixedParams, fixedParamValues,
00103 tuple(fixedFields), tuple(fixedFieldValues)));
00104
00105 assembler_ = rcp(new Assembler(mesh, eqnSet, tuple(vecType), tuple(vecType), partitionBCs));
00106
00107 VectorSpace<double> domain = assembler_->solnVecSpace();
00108 VectorSpace<double> range = assembler_->rowVecSpace();
00109
00110 setDomainAndRange(domain, range);
00111 }
00112
00113
00114
00115 NLOp::NLOp(const Mesh& mesh,
00116 const Expr& eqn,
00117 const Expr& bc,
00118 const BlockArray& test,
00119 const BlockArray& unk,
00120 const Expr& u0)
00121 : NonlinearOperatorBase<double>(),
00122 assembler_(),
00123 u0_(u0),
00124 params_()
00125 {
00126 TimeMonitor timer(nlpCtorTimer());
00127 bool partitionBCs = false;
00128
00129 Array<Expr> v(test.size());
00130 Array<Expr> u(unk.size());
00131 Array<Expr> u0Array(unk.size());
00132
00133 Array<VectorType<double> > testVecType(test.size());
00134 Array<VectorType<double> > unkVecType(unk.size());
00135
00136 for (int i=0; i<test.size(); i++)
00137 {
00138 v[i] = test[i].expr().flattenSpectral();
00139 testVecType[i] = test[i].vecType();
00140 }
00141
00142 for (int i=0; i<unk.size(); i++)
00143 {
00144 u[i] = unk[i].expr().flattenSpectral();
00145 u0Array[i] = u0[i].flattenSpectral();
00146 unkVecType[i] = unk[i].vecType();
00147 }
00148
00149 Expr unkParams;
00150 Expr fixedParams;
00151 Expr fixedFields;
00152 Expr unkParamValues;
00153 Expr fixedParamValues;
00154 Expr fixedFieldValues;
00155
00156 RCP<EquationSet> eqnSet
00157 = rcp(new EquationSet(eqn, bc,
00158 v, u, u0Array,
00159 unkParams, unkParamValues,
00160 fixedParams, fixedParamValues,
00161 tuple(fixedFields), tuple(fixedFieldValues)));
00162
00163 assembler_ = rcp(new Assembler(mesh, eqnSet, testVecType, unkVecType, partitionBCs));
00164
00165 VectorSpace<double> domain = assembler_->solnVecSpace();
00166 VectorSpace<double> range = assembler_->rowVecSpace();
00167
00168 setDomainAndRange(domain, range);
00169 }
00170
00171
00172
00173 NLOp::NLOp(const Mesh& mesh,
00174 const Expr& eqn,
00175 const Expr& bc,
00176 const Expr& test,
00177 const Expr& unk,
00178 const Expr& u0,
00179 const Expr& params,
00180 const Expr& paramValues,
00181 const VectorType<double>& vecType,
00182 bool partitionBCs)
00183 : NonlinearOperatorBase<double>(),
00184 assembler_(),
00185 J_(),
00186 u0_(u0),
00187 params_(params)
00188 {
00189 TimeMonitor timer(nlpCtorTimer());
00190
00191 Expr unkParams;
00192 Expr fixedFields;
00193 Expr unkParamValues;
00194 Expr fixedFieldValues;
00195
00196 RCP<EquationSet> eqnSet
00197 = rcp(new EquationSet(
00198 eqn, bc, tuple(test.flattenSpectral()),
00199 tuple(unk.flattenSpectral()), tuple(u0),
00200 unkParams, unkParamValues,
00201 params, paramValues,
00202 tuple(fixedFields), tuple(fixedFieldValues)));
00203
00204 assembler_ = rcp(new Assembler(mesh, eqnSet, tuple(vecType), tuple(vecType), partitionBCs));
00205
00206 VectorSpace<double> domain = assembler_->solnVecSpace();
00207 VectorSpace<double> range = assembler_->rowVecSpace();
00208
00209 setDomainAndRange(domain, range);
00210 }
00211
00212
00213 NLOp::NLOp(const RCP<Assembler>& assembler,
00214 const Expr& u0)
00215 : NonlinearOperatorBase<double>(),
00216 assembler_(assembler),
00217 J_(),
00218 u0_(u0),
00219 params_()
00220 {
00221 TimeMonitor timer(nlpCtorTimer());
00222
00223 VectorSpace<double> domain = assembler_->solnVecSpace();
00224 VectorSpace<double> range = assembler_->rowVecSpace();
00225
00226 setDomainAndRange(domain, range);
00227 }
00228
00229
00230 void NLOp::updateDiscreteFunctionValue(const Vector<double>& vec) const
00231 {
00232 setDiscreteFunctionVector(u0_, vec.copy());
00233 }
00234
00235 Vector<double> NLOp::getInitialGuess() const
00236 {
00237 return getDiscreteFunctionVector(u0_);
00238 }
00239
00240 void NLOp::setInitialGuess(const Expr& u0New)
00241 {
00242 Vector<double> vec = getDiscreteFunctionVector(u0New);
00243 setEvalPt(vec);
00244 updateDiscreteFunctionValue(vec);
00245 }
00246
00247
00248 LinearOperator<double> NLOp::allocateJacobian() const
00249 {
00250 return assembler_->allocateMatrix();
00251 }
00252
00253
00254 LinearOperator<double> NLOp::
00255 computeJacobianAndFunction(Vector<double>& functionValue) const
00256 {
00257 updateDiscreteFunctionValue(currentEvalPt());
00258
00259 Array<Vector<double> > mv(1);
00260 mv[0].acceptCopyOf(functionValue);
00261 assembler_->assemble(J_, mv);
00262 functionValue.acceptCopyOf(mv[0]);
00263
00264 return J_;
00265 }
00266
00267 void NLOp::
00268 computeJacobianAndFunction(LinearOperator<double>& J,
00269 Vector<double>& resid) const
00270 {
00271 TEUCHOS_TEST_FOR_EXCEPTION(currentEvalPt().ptr().get()==0, std::runtime_error,
00272 "null evaluation point in "
00273 "NLOp::jacobian()");
00274
00275 TEUCHOS_TEST_FOR_EXCEPTION(J.ptr().get()==0, std::runtime_error,
00276 "null Jacobian pointer in "
00277 "NLOp::jacobian()");
00278
00279 TEUCHOS_TEST_FOR_EXCEPTION(resid.ptr().get()==0, std::runtime_error,
00280 "null residual pointer in "
00281 "NLOp::jacobian()");
00282
00283 Array<Vector<double> > mv(1);
00284 mv[0] = resid;
00285
00286 updateDiscreteFunctionValue(currentEvalPt());
00287 assembler_->assemble(J, mv);
00288
00289 resid.acceptCopyOf(mv[0]);
00290
00291 J_ = J;
00292 }
00293
00294
00295
00296
00297 Vector<double> NLOp::computeFunctionValue() const
00298 {
00299 TEUCHOS_TEST_FOR_EXCEPTION(currentEvalPt().ptr().get()==0, std::runtime_error,
00300 "null evaluation point in "
00301 "NLOp::computeFunctionValue()");
00302
00303 VectorSpace<double> spc = range();
00304 Vector<double> rtn = spc.createMember();
00305
00306 Array<Vector<double> > mv(1);
00307 mv[0] = rtn;
00308
00309 updateDiscreteFunctionValue(currentEvalPt());
00310 assembler_->assemble(mv);
00311
00312 rtn.acceptCopyOf(mv[0]);
00313
00314 return rtn;
00315 }
00316
00317
00318
00319 void NLOp::computeFunctionValue(Vector<double>& resid) const
00320 {
00321
00322
00323
00324 TEUCHOS_TEST_FOR_EXCEPTION(currentEvalPt().ptr().get()==0, std::runtime_error,
00325 "null evaluation point in "
00326 "NLOp::computeFunctionValue()");
00327
00328 TEUCHOS_TEST_FOR_EXCEPTION(resid.ptr().get()==0, std::runtime_error,
00329 "null residual pointer in "
00330 "NLOp::computeFunctionValue()");
00331
00332 updateDiscreteFunctionValue(currentEvalPt());
00333
00334
00335 Array<Vector<double> > mv(1);
00336 mv[0] = resid;
00337
00338 assembler_->assemble(mv);
00339
00340 resid.acceptCopyOf(mv[0]);
00341 }
00342
00343
00344
00345 Expr NLOp::
00346 computeSensitivities(const LinearSolver<double>& solver) const
00347 {
00348 TEUCHOS_TEST_FOR_EXCEPTION(currentEvalPt().ptr().get()==0, std::runtime_error,
00349 "null evaluation point in "
00350 "NLOp::computeSensitivities()");
00351
00352 TEUCHOS_TEST_FOR_EXCEPTION(J_.ptr().get()==0, std::runtime_error,
00353 "null Jacobian pointer in "
00354 "NLOp::computeSensitivities()");
00355
00356 TEUCHOS_TEST_FOR_EXCEPTION(params_.ptr().get()==0, std::runtime_error,
00357 "null parameters in NLOp::computeSensitivities()");
00358
00359 updateDiscreteFunctionValue(currentEvalPt());
00360
00361 Array<Vector<double> > mv(params_.size());
00362
00363 assembler_->assembleSensitivities(J_, mv);
00364
00365 LinearSolveDriver driver;
00366
00367 Expr sens;
00368 int vrb = 0;
00369 Array<Array<string> > names(params_.size());
00370 for (int i=0; i<params_.size(); i++)
00371 {
00372 names[i].resize(u0_.size());
00373 for (int j=0; j<u0_.size(); j++)
00374 names[i][j]="sens(" + u0_[j].toString() + ", " + params_[i].toString() + ")";
00375 mv[i].scale(-1.0);
00376 }
00377
00378 driver.solve(solver, J_, mv, assembler_->solutionSpace(), names, vrb, sens);
00379 return sens;
00380 }
00381
00382
00383 void NLOp::reAssembleProblem() const
00384 {
00385 assembler_->flushConfiguration();
00386 }