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