00001 #include "PDEOptPDEConstrainedObjBase.hpp"
00002 #include "Sundance.hpp"
00003
00004
00005 using namespace Teuchos;
00006 using namespace Playa;
00007 using namespace Sundance;
00008
00009 using std::endl;
00010 using std::runtime_error;
00011
00012
00013 PDEConstrainedObjBase::PDEConstrainedObjBase(
00014 const Functional& lagrangian,
00015 const Array<Expr>& stateVarVals,
00016 const Array<Expr>& adjointVarVals,
00017 const Expr& designVarVal,
00018 int verb)
00019 : ObjectiveBase(verb),
00020 Lagrangian_(lagrangian),
00021 designVarVal_(designVarVal),
00022 stateVarVals_(stateVarVals),
00023 adjointVarVals_(adjointVarVals),
00024 fEval_(),
00025 numFuncEvals_(0),
00026 numGradEvals_(0),
00027 invHScale_(1.0),
00028 iterCallback_()
00029 {}
00030
00031
00032 PDEConstrainedObjBase::PDEConstrainedObjBase(
00033 const Functional& lagrangian,
00034 const Array<Expr>& stateVarVals,
00035 const Array<Expr>& adjointVarVals,
00036 const Expr& designVarVal,
00037 const RCP<IterCallbackBase>& iterCallback,
00038 int verb)
00039 : ObjectiveBase(verb),
00040 Lagrangian_(lagrangian),
00041 designVarVal_(designVarVal),
00042 stateVarVals_(stateVarVals),
00043 adjointVarVals_(adjointVarVals),
00044 fEval_(),
00045 numFuncEvals_(0),
00046 numGradEvals_(0),
00047 invHScale_(1.0),
00048 iterCallback_(iterCallback)
00049 {}
00050
00051
00052 void PDEConstrainedObjBase::init(
00053 const Array<Expr>& stateVars,
00054 const Array<Expr>& adjointVars,
00055 const Expr& designVar)
00056 {
00057 TEUCHOS_TEST_FOR_EXCEPTION(stateVars.size() != adjointVars.size(),
00058 RuntimeError,
00059 "number of state and adjoint variables should be identical");
00060
00061 TEUCHOS_TEST_FOR_EXCEPTION(stateVars.size() != stateVarVals_.size(),
00062 RuntimeError,
00063 "number of state variables and state values "
00064 "should be identical");
00065
00066 TEUCHOS_TEST_FOR_EXCEPTION(adjointVars.size() != adjointVarVals_.size(),
00067 RuntimeError,
00068 "number of adjoint variables and adjoint values "
00069 "should be identical");
00070
00071 Array<Array<Expr> > fixedVarsInStateEqns(stateVars.size());
00072 Array<Array<Expr> > fixedVarsInStateEqnsVals(stateVars.size());
00073
00074 Array<Array<Expr> > fixedVarsInAdjointEqns(stateVars.size());
00075 Array<Array<Expr> > fixedVarsInAdjointEqnsVals(stateVars.size());
00076
00077 for (int i=0; i<stateVars.size(); i++)
00078 {
00079
00080 fixedVarsInStateEqns[i].append(designVar);
00081 fixedVarsInStateEqnsVals[i].append(designVarVal_);
00082 for (int j=0; j<stateVars.size(); j++)
00083 {
00084 if (i != j)
00085 {
00086 fixedVarsInStateEqns[i].append(stateVars[j]);
00087 fixedVarsInStateEqnsVals[i].append(stateVarVals_[j]);
00088 fixedVarsInAdjointEqns[i].append(adjointVars[j]);
00089 fixedVarsInAdjointEqnsVals[i].append(adjointVarVals_[j]);
00090 fixedVarsInStateEqns[i].append(adjointVars[j]);
00091 fixedVarsInStateEqnsVals[i].append(adjointVarVals_[j]);
00092 fixedVarsInAdjointEqns[i].append(stateVars[j]);
00093 fixedVarsInAdjointEqnsVals[i].append(stateVarVals_[j]);
00094 }
00095 }
00096 }
00097
00098 initEquations(stateVars, adjointVars,
00099 fixedVarsInStateEqns, fixedVarsInStateEqnsVals,
00100 fixedVarsInAdjointEqns, fixedVarsInAdjointEqnsVals);
00101
00102 Array<Expr> allVarArray = stateVars;
00103 Array<Expr> allVarValArray = stateVarVals_;
00104 for (int i=0; i<adjointVars.size(); i++)
00105 {
00106 allVarArray.append(adjointVars[i]);
00107 allVarValArray.append(adjointVarVals_[i]);
00108 }
00109 Expr allVars = new ListExpr(allVarArray);
00110 Expr allVarVals = new ListExpr(allVarValArray);
00111
00112 fEval_ = Lagrangian_.evaluator(designVar, designVarVal_,
00113 allVars, allVarVals);
00114 }
00115
00116
00117 void PDEConstrainedObjBase::evalGrad(const Vector<double>& x, double& f,
00118 Vector<double>& grad) const
00119 {
00120 Tabs tabs(0);
00121 PLAYA_MSG2(verb(), tabs << "in evalGrad()");
00122
00123
00124 solveStateAndAdjoint(x);
00125
00126
00127 Expr gradExpr = fEval_.evalGradient(f);
00128 grad = getDiscreteFunctionVector(gradExpr).copy();
00129
00130 PLAYA_MSG2(verb(), tabs << "function value = " << f);
00131 PLAYA_MSG5(verb(), tabs << "|gradient| is " << grad.norm2());
00132 PLAYA_MSG5(verb(), tabs << "gradient is " << grad);
00133 numFuncEvals_++;
00134 numGradEvals_++;
00135 }
00136
00137 void PDEConstrainedObjBase::eval(const Vector<double>& x, double& f) const
00138 {
00139 Tabs tabs(0);
00140 PLAYA_MSG2(verb(), tabs << "in eval()");
00141
00142
00143 solveState(x);
00144
00145
00146
00147
00148 for (int i=0; i<adjointVarVals_.size(); i++)
00149 {
00150
00151 Vector<double> adjVec = getDiscreteFunctionVector(adjointVarVals_[i]);
00152 adjVec.zero();
00153 setDiscreteFunctionVector(adjointVarVals_[i], adjVec);
00154 }
00155
00156
00157 f = fEval_.evaluate();
00158 PLAYA_MSG2(verb(), tabs << "function value = " << f);
00159 numFuncEvals_++;
00160 }
00161
00162
00163 Vector<double> PDEConstrainedObjBase::getInit() const
00164 {
00165 return getDiscreteFunctionVector(designVarVal());
00166 }
00167
00168 void PDEConstrainedObjBase:: iterationCallback(const Vector<double>& x,
00169 int iter) const
00170 {
00171 if (iterCallback_.ptr().get() != 0)
00172 {
00173 iterCallback_->call(this, iter);
00174 }
00175 }
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186