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 "SundanceFunctionalEvaluator.hpp"
00032 #include "SundanceEquationSet.hpp"
00033 #include "SundanceAssembler.hpp"
00034 #include "SundanceSymbPreprocessor.hpp"
00035 #include "SundanceDiscreteFunction.hpp"
00036 #include "SundanceEquationSet.hpp"
00037 #include "SundanceDiscreteSpace.hpp"
00038 #include "SundanceOut.hpp"
00039 #include "PlayaTabs.hpp"
00040
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043 using namespace Playa;
00044
00045
00046 using std::endl;
00047 using std::setw;
00048
00049
00050 namespace Sundance
00051 {
00052 double evaluateIntegral(const Mesh& mesh, const Expr& expr)
00053 {
00054 FunctionalEvaluator eval(mesh, expr);
00055 return eval.evaluate();
00056 }
00057 }
00058
00059 FunctionalEvaluator::FunctionalEvaluator()
00060 : assembler_(),
00061 varValues_(),
00062 vecType_(),
00063 gradient_()
00064 {}
00065
00066 FunctionalEvaluator::FunctionalEvaluator(const Mesh& mesh,
00067 const Expr& integral)
00068 : assembler_(),
00069 varValues_(),
00070 vecType_(),
00071 gradient_(1)
00072 {
00073 Array<Expr> fields;
00074 Expr bcs;
00075 Expr params;
00076
00077
00078 RCP<EquationSet> eqnSet
00079 = rcp(new EquationSet(integral, bcs, params, params, fields, fields));
00080
00081
00082 assembler_ = rcp(new Assembler(mesh, eqnSet));
00083 }
00084
00085
00086 FunctionalEvaluator::FunctionalEvaluator(const Mesh& mesh,
00087 const Expr& integral,
00088 const Expr& bcs,
00089 const Expr& var,
00090 const Expr& varValues,
00091 const VectorType<double>& vectorType)
00092 : assembler_(),
00093 varValues_(varValues),
00094 vecType_(vectorType),
00095 gradient_(1)
00096 {
00097 Array<Expr> v = tuple(var.flatten());
00098 Array<Expr> v0 = tuple(varValues.flatten());
00099 Array<Expr> fixed;
00100 Expr params;
00101
00102 RCP<EquationSet> eqnSet
00103 = rcp(new EquationSet(integral, bcs, v, v0, params, params, fixed, fixed));
00104
00105 assembler_ = rcp(new Assembler(mesh, eqnSet, tuple(vectorType), tuple(vectorType), false));
00106 }
00107
00108
00109 FunctionalEvaluator::FunctionalEvaluator(const Mesh& mesh,
00110 const Expr& integral,
00111 const Expr& bcs,
00112 const Expr& vars,
00113 const Expr& varEvalPts,
00114 const Expr& fields,
00115 const Expr& fieldValues,
00116 const VectorType<double>& vectorType)
00117 : assembler_(),
00118 varValues_(varEvalPts),
00119 vecType_(vectorType),
00120 gradient_(1)
00121 {
00122 Array<Expr> f = tuple(fields.flatten());
00123 Array<Expr> f0 = tuple(fieldValues.flatten());
00124 Array<Expr> v = tuple(vars.flatten());
00125 Array<Expr> v0 = tuple(varEvalPts.flatten());
00126
00127 Expr params;
00128
00129 RCP<EquationSet> eqnSet
00130 = rcp(new EquationSet(integral, bcs, v, v0, params, params, f, f0));
00131
00132 assembler_ = rcp(new Assembler(mesh, eqnSet, tuple(vectorType), tuple(vectorType), false));
00133 }
00134
00135 double FunctionalEvaluator::evaluate() const
00136 {
00137 double value;
00138 assembler_->evaluate(value);
00139 return value;
00140 }
00141
00142
00143 Vector<double> FunctionalEvaluator::evalGradientVector(double& value) const
00144 {
00145 assembler_->evaluate(value, gradient_);
00146
00147 return gradient_[0];
00148 }
00149
00150 Expr FunctionalEvaluator::evalGradient(double& value) const
00151 {
00152 Vector<double> g = evalGradientVector(value);
00153
00154 Array<Expr> rtn(assembler_->rowSpace().size());
00155 for (int i=0; i<rtn.size(); i++)
00156 {
00157 std::string name = "gradient";
00158 if (rtn.size() > 1) name += "[" + Teuchos::toString(i) + "]";
00159 rtn[i] = new DiscreteFunction(*(assembler_->rowSpace()[i]),
00160 g.getBlock(i), name);
00161 }
00162 if ((int)rtn.size()==1)
00163 {
00164 return rtn[0];
00165 }
00166 else
00167 {
00168 return new ListExpr(rtn);
00169 }
00170 }
00171
00172
00173 double FunctionalEvaluator::fdGradientCheck(double h) const
00174 {
00175 bool showAll = false;
00176 Tabs tabs;
00177 double f0, fPlus, fMinus;
00178 Expr gradF0 = evalGradient(f0);
00179
00180 FancyOStream& os = Out::os();
00181
00182
00183 DiscreteFunction* df = DiscreteFunction::discFunc(varValues_);
00184 DiscreteFunction* dg = DiscreteFunction::discFunc(gradF0);
00185 Vector<double> x = df->getVector();
00186 Vector<double> x0 = x.copy();
00187 Vector<double> gf = dg->getVector();
00188
00189 RCP<GhostView<double> > xView = df->ghostView();
00190 RCP<GhostView<double> > gradF0View = dg->ghostView();
00191
00192
00193 TEUCHOS_TEST_FOR_EXCEPTION(xView.get() == 0, std::runtime_error,
00194 "bad pointer in FunctionalEvaluator::fdGradientCheck");
00195 TEUCHOS_TEST_FOR_EXCEPTION(gradF0View.get() == 0, std::runtime_error,
00196 "bad pointer in FunctionalEvaluator::fdGradientCheck");
00197
00198 int nTot = x.space().dim();
00199 int n = x.space().numLocalElements();
00200 int lowestIndex = x.space().baseGlobalNaturalIndex();
00201
00202 os << tabs << "doing fd check: h=" << h << std::endl;
00203 Array<double> df_dx(n);
00204
00205 int localIndex = 0;
00206 for (int globalIndex=0; globalIndex<nTot; globalIndex++)
00207 {
00208 double tmp=0.0;
00209 bool isLocal = globalIndex >= lowestIndex
00210 && globalIndex < (lowestIndex+n);
00211 if (isLocal)
00212 {
00213 tmp = xView->getElement(globalIndex);
00214 loadable(x)->setElement(globalIndex, tmp + h);
00215 }
00216
00217 df->setVector(x);
00218 fPlus = evaluate();
00219 if (isLocal)
00220 {
00221 loadable(x)->setElement(globalIndex, tmp - h);
00222 }
00223
00224 df->setVector(x);
00225 fMinus = evaluate();
00226
00227 if (isLocal)
00228 {
00229 df_dx[localIndex] = (fPlus - fMinus)/2.0/h;
00230 os << "g=" << setw(5) << globalIndex << ", l=" << setw(5) << localIndex << " f0="
00231 << setw(12) << f0
00232 << " fPlus=" << setw(12) << fPlus << " fMinus=" << setw(12) << fMinus << " df_dx="
00233 << setw(12) << df_dx[localIndex] << std::endl;
00234 if (showAll)
00235 {
00236 os << "i " << globalIndex << " x_i=" << tmp
00237 << " f(x)=" << f0
00238 << " f(x+h)=" << fPlus
00239 << " f(x-h)=" << fMinus << std::endl;
00240 }
00241 loadable(x)->setElement(globalIndex, tmp);
00242 localIndex++;
00243 }
00244 df->setVector(x);
00245 }
00246
00247 double localMaxErr = 0.0;
00248
00249 showAll = true;
00250 VectorSpace<double> space = x.space();
00251 for (int i=0; i<space.numLocalElements(); i++)
00252 {
00253 double num = fabs(df_dx[i]-gf[i]);
00254 double den = fabs(df_dx[i]) + fabs(gf[i]) + 1.0e-14;
00255 double r = 0.0;
00256 if (fabs(den) > 1.0e-16) r = num/den;
00257 else r = 1.0;
00258 if (showAll)
00259 {
00260 os << "i " << i;
00261 os << " FD=" << df_dx[i]
00262 << " grad=" << gf[i]
00263 << " r=" << r << std::endl;
00264 }
00265 if (localMaxErr < r) localMaxErr = r;
00266 }
00267 os << "local max error = " << localMaxErr << std::endl;
00268
00269 double maxErr = localMaxErr;
00270 df->mesh().comm().allReduce((void*) &localMaxErr, (void*) &maxErr, 1,
00271 MPIDataType::doubleType(), MPIOp::maxOp());
00272 os << tabs << "fd check: max error = " << maxErr << std::endl;
00273
00274 return maxErr;
00275 }
00276