00001 #include "PDEOptNonlinearPDEConstrainedObj.hpp"
00002 #include "Sundance.hpp"
00003
00004 using namespace Teuchos;
00005 using namespace Playa;
00006
00007 namespace Sundance
00008 {
00009
00010 NonlinearPDEConstrainedObj::NonlinearPDEConstrainedObj(
00011 const Functional& lagrangian,
00012 const Expr& stateVars,
00013 const Expr& stateVarVals,
00014 const Expr& adjointVars,
00015 const Expr& adjointVarVals,
00016 const Expr& designVar,
00017 const Expr& designVarVal,
00018 const NOXSolver& solver,
00019 const LinearSolver<double>& adjSolver,
00020 int verb)
00021 : PDEConstrainedObjBase(lagrangian, tuple(stateVarVals),
00022 tuple(adjointVarVals), designVarVal, verb),
00023 stateProbs_(),
00024 adjointProbs_(),
00025 solver_(solver),
00026 adjSolver_(adjSolver)
00027 {
00028 init(tuple(stateVars), tuple(adjointVars), designVar);
00029 }
00030
00031
00032 NonlinearPDEConstrainedObj::NonlinearPDEConstrainedObj(
00033 const Functional& lagrangian,
00034 const Array<Expr>& stateVars,
00035 const Array<Expr>& stateVarVals,
00036 const Array<Expr>& adjointVars,
00037 const Array<Expr>& adjointVarVals,
00038 const Expr& designVar,
00039 const Expr& designVarVal,
00040 const NOXSolver& solver,
00041 const LinearSolver<double>& adjSolver,
00042 int verb)
00043 : PDEConstrainedObjBase(lagrangian, stateVarVals,
00044 adjointVarVals, designVarVal, verb),
00045 stateProbs_(),
00046 adjointProbs_(),
00047 solver_(solver),
00048 adjSolver_(adjSolver)
00049 {
00050 init(stateVars, adjointVars, designVar);
00051 }
00052
00053 NonlinearPDEConstrainedObj::NonlinearPDEConstrainedObj(
00054 const Functional& lagrangian,
00055 const Expr& stateVars,
00056 const Expr& stateVarVals,
00057 const Expr& adjointVars,
00058 const Expr& adjointVarVals,
00059 const Expr& designVar,
00060 const Expr& designVarVal,
00061 const NOXSolver& solver,
00062 const LinearSolver<double>& adjSolver,
00063 const RCP<IterCallbackBase>& iterCallback,
00064 int verb)
00065 : PDEConstrainedObjBase(lagrangian, tuple(stateVarVals),
00066 tuple(adjointVarVals), designVarVal, iterCallback, verb),
00067 stateProbs_(),
00068 adjointProbs_(),
00069 solver_(solver),
00070 adjSolver_(adjSolver)
00071 {
00072 init(tuple(stateVars), tuple(adjointVars), designVar);
00073 }
00074
00075
00076 NonlinearPDEConstrainedObj::NonlinearPDEConstrainedObj(
00077 const Functional& lagrangian,
00078 const Array<Expr>& stateVars,
00079 const Array<Expr>& stateVarVals,
00080 const Array<Expr>& adjointVars,
00081 const Array<Expr>& adjointVarVals,
00082 const Expr& designVar,
00083 const Expr& designVarVal,
00084 const NOXSolver& solver,
00085 const LinearSolver<double>& adjSolver,
00086 const RCP<IterCallbackBase>& iterCallback,
00087 int verb)
00088 : PDEConstrainedObjBase(lagrangian, stateVarVals,
00089 adjointVarVals, designVarVal, iterCallback, verb),
00090 stateProbs_(),
00091 adjointProbs_(),
00092 solver_(solver),
00093 adjSolver_(adjSolver)
00094 {
00095 init(stateVars, adjointVars, designVar);
00096 }
00097
00098
00099 void NonlinearPDEConstrainedObj::initEquations(
00100 const Array<Expr>& stateVars,
00101 const Array<Expr>& adjointVars,
00102 const Array<Array<Expr> >& fixedVarsInStateEqns,
00103 const Array<Array<Expr> >& fixedVarsInStateEqnsVals,
00104 const Array<Array<Expr> >& fixedVarsInAdjointEqns,
00105 const Array<Array<Expr> >& fixedVarsInAdjointEqnsVals
00106 )
00107 {
00108 Tabs tab(0);
00109 PLAYA_MSG2(verb(), tab << "setting up nonlinear equations");
00110
00111 for (int i=0; i<stateVars.size(); i++)
00112 {
00113 Tabs tab1;
00114 PLAYA_MSG3(verb(), tab1 << "setting up state equation #" << i);
00115 Expr fixedVars = new ListExpr(fixedVarsInStateEqns[i]);
00116 Expr fixedVarVals = new ListExpr(fixedVarsInStateEqnsVals[i]);
00117 PLAYA_MSG3(verb(), tab1 << "Fixed vars are: " << fixedVars);
00118 NonlinearProblem stateProb
00119 = Lagrangian().nonlinearVariationalProb(adjointVars[i],
00120 adjointVarVals(i),
00121 stateVars[i], stateVarVals(i),
00122 fixedVars, fixedVarVals);
00123
00124 stateProbs_.append(stateProb);
00125 }
00126
00127 for (int i=0; i<adjointVars.size(); i++)
00128 {
00129 Tabs tab1;
00130 PLAYA_MSG3(verb(), tab1 << "setting up adjoint equation #" << i);
00131 Expr fixedVars = new ListExpr(fixedVarsInAdjointEqns[i]);
00132 Expr fixedVarVals = new ListExpr(fixedVarsInAdjointEqnsVals[i]);
00133 PLAYA_MSG3(verb(), tab1 << "Fixed vars are: " << fixedVars);
00134 LinearProblem adjointProb
00135 = Lagrangian().linearVariationalProb(stateVars[i], stateVarVals(i),
00136 adjointVars[i],
00137 fixedVars, fixedVarVals);
00138
00139 adjointProbs_.append(adjointProb);
00140 }
00141
00142 PLAYA_MSG2(verb(), tab << "done setting up nonlinear equations");
00143 }
00144
00145
00146
00147
00148 void NonlinearPDEConstrainedObj::solveState(const Vector<double>& x) const
00149 {
00150 Tabs tab(0);
00151 PLAYA_MSG2(verb(), tab << "solving state");
00152 PLAYA_MSG3(verb(), tab << "|x|=" << x.norm2());
00153 PLAYA_MSG5(verb(), tab << "x=" << endl << tab << x.norm2());
00154 setDiscreteFunctionVector(designVarVal(), x);
00155
00156
00157 for (int i=0; i<stateProbs_.size(); i++)
00158 {
00159 PLAYA_MSG3(verb(), tab << "state eqn=" << i);
00160 SolverState<double> status
00161 = stateProbs_[i].solve(solver_);
00162 TEUCHOS_TEST_FOR_EXCEPTION(status.finalState() != SolveConverged,
00163 std::runtime_error,
00164 "state equation could not be solved: status="
00165 << status);
00166 }
00167
00168 PLAYA_MSG2(verb(), tab << "done state solve");
00169
00170 statePostprocCallback();
00171 }
00172
00173
00174
00175 void NonlinearPDEConstrainedObj
00176 ::solveStateAndAdjoint(const Vector<double>& x) const
00177 {
00178 Tabs tab(0);
00179 PLAYA_MSG2(verb(), tab << "solving state and adjoint");
00180 PLAYA_MSG3(verb(), tab << "|x|=" << x.norm2());
00181 PLAYA_MSG5(verb(), tab << "x=" << endl << tab << x.norm2());
00182
00183 Tabs tab1;
00184 setDiscreteFunctionVector(designVarVal(), x);
00185
00186 PLAYA_MSG3(verb(), tab1 << "solving state eqns");
00187
00188 for (int i=0; i<stateProbs_.size(); i++)
00189 {
00190 PLAYA_MSG3(verb(), tab << "state eqn=" << i);
00191 SolverState<double> status
00192 = stateProbs_[i].solve(solver_);
00193
00194
00195
00196 bool dumpBadSolve = false;
00197 if (dumpBadSolve && status.finalState() != SolveConverged)
00198 {
00199 FieldWriter w = new VTKWriter("badSolve");
00200 w.addMesh(Lagrangian().mesh());
00201 w.addField("designVar", new ExprFieldWrapper(designVarVal()));
00202 for (int j=0; j<i; j++)
00203 {
00204 Expr tmp = stateVarVals(j).flatten();
00205 for (int k=0; k<tmp.size(); k++)
00206 {
00207 w.addField("stateVar-"+Teuchos::toString(j)+"-"+Teuchos::toString(k),
00208 new ExprFieldWrapper(tmp[k]));
00209 }
00210 }
00211 w.write();
00212 }
00213 TEUCHOS_TEST_FOR_EXCEPTION(status.finalState() != SolveConverged,
00214 std::runtime_error,
00215 "state equation " << i
00216 << " could not be solved: status="
00217 << status);
00218 }
00219
00220 PLAYA_MSG3(verb(), tab1 << "done solving state eqns");
00221
00222
00223 statePostprocCallback();
00224
00225 PLAYA_MSG3(verb(), tab1 << "solving adjoint eqns");
00226
00227
00228 for (int i=adjointProbs_.size()-1; i>=0; i--)
00229 {
00230 PLAYA_MSG3(verb(), tab << "adjoint eqn=" << i);
00231 SolverState<double> status
00232 = adjointProbs_[i].solve(adjSolver_, adjointVarVals(i));
00233
00234
00235
00236 bool dumpBadSolve = false;
00237 if (dumpBadSolve && status.finalState() != SolveConverged)
00238 {
00239 FieldWriter w = new VTKWriter("badSolve");
00240 w.addMesh(Lagrangian().mesh());
00241 w.addField("designVar", new ExprFieldWrapper(designVarVal()));
00242 for (int j=0; j<stateProbs_.size(); j++)
00243 {
00244 Expr tmp = stateVarVals(j).flatten();
00245 for (int k=0; k<tmp.size(); k++)
00246 {
00247 w.addField("stateVar-"+Teuchos::toString(j)+"-"+Teuchos::toString(k),
00248 new ExprFieldWrapper(tmp[k]));
00249 }
00250 }
00251 for (int j=adjointProbs_.size()-1; j>i; j--)
00252 {
00253 Expr tmp = adjointVarVals(j).flatten();
00254 for (int k=0; k<tmp.size(); k++)
00255 {
00256 w.addField("adjointVar-"+Teuchos::toString(j)+"-"+Teuchos::toString(k),
00257 new ExprFieldWrapper(tmp[k]));
00258 }
00259
00260 }
00261 w.write();
00262
00263 }
00264 TEUCHOS_TEST_FOR_EXCEPTION(status.finalState() != SolveConverged,
00265 std::runtime_error,
00266 "adjoint equation " << i
00267 << " could not be solved: status="
00268 << status.stateDescription());
00269 }
00270 PLAYA_MSG3(verb(), tab1 << "done solving adjoint eqns");
00271 PLAYA_MSG2(verb(), tab1 << "done solving state and adjoint eqns");
00272 }
00273
00274
00275
00276
00277 }
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
