00001
00002
00003
00004 #include "Sundance.hpp"
00005
00006
00007
00008
00009
00010
00011
00012
00013 const double pi = 4.0*atan(1.0);
00014
00015
00016
00017 class WestEdgeTest : public CellPredicateFunctorBase,
00018 public Playa::Handleable<CellPredicateFunctorBase>
00019 {
00020 public:
00021 WestEdgeTest() : CellPredicateFunctorBase("WestEdgeTest") {}
00022 virtual ~WestEdgeTest() {}
00023 virtual bool operator()(const Point& x) const {return fabs(x[0])<1.0e-10;}
00024 GET_RCP(CellPredicateFunctorBase);
00025 };
00026
00027
00028 CELL_PREDICATE(EastEdgeTest, {return fabs(x[0]-1.0) < 1.0e-10;})
00029 CELL_PREDICATE(SouthEdgeTest, {return fabs(x[1]) < 1.0e-10;})
00030 CELL_PREDICATE(NorthEdgeTest, {return fabs(x[1]-1.0) < 1.0e-10;})
00031
00032 bool inlineNewtonSolve(NonlinearProblem prob,
00033 Expr uNewt,
00034 int maxIters,
00035 double newtTol);
00036
00037 bool noxNewtonSolve(const NonlinearProblem& prob, Expr uNewt);
00038
00039 bool playaNewtonArmijoSolve(const NonlinearProblem& prob, Expr uNewt);
00040
00041 int main(int argc, char** argv)
00042 {
00043 try
00044 {
00045 Sundance::init(&argc, &argv);
00046
00047
00048 VectorType<double> vecType = new EpetraVectorType();
00049
00050
00051
00052 MeshType meshType = new BasicSimplicialMeshType();
00053 int nx = 128;
00054 MeshSource mesher = new PartitionedRectangleMesher(
00055 0.0, 1.0, nx, 1,
00056 0.0, 1.0, nx, 1, meshType);
00057 Mesh mesh = mesher.getMesh();
00058
00059
00060
00061 CellFilter interior = new MaximalCellFilter();
00062 CellFilter edges = new DimensionalCellFilter(1);
00063 CellFilter northEdge = edges.subset(new NorthEdgeTest());
00064 CellFilter southEdge = edges.subset(new SouthEdgeTest());
00065 CellFilter eastEdge = edges.subset(new EastEdgeTest());
00066 CellFilter westEdge = edges.subset(new WestEdgeTest());
00067
00068
00069
00070
00071 BasisFamily bas = new Lagrange(1);
00072 Expr u = new UnknownFunction(bas, "u");
00073 Expr v = new TestFunction(bas, "v");
00074
00075
00076 Expr dx = new Derivative(0);
00077 Expr dy = new Derivative(1);
00078 Expr grad = List(dx, dy);
00079 Expr x = new CoordExpr(0);
00080 Expr y = new CoordExpr(1);
00081
00082
00083 DiscreteSpace discSpace(mesh, bas, vecType);
00084 Expr uNewt = new DiscreteFunction(discSpace, 1.0, "uNewt");
00085
00086
00087 QuadratureFamily quad = new GaussianQuadrature(4);
00088 WatchFlag watchMe("watch me");
00089 watchMe.setParam("symbolic preprocessing", 1);
00090 watchMe.setParam("discrete function evaluation", 3);
00091 watchMe.setParam("integration setup", 6);
00092 watchMe.setParam("integration", 6);
00093 watchMe.setParam("fill", 6);
00094 watchMe.setParam("evaluation", 2);
00095 watchMe.deactivate();
00096
00097 Expr eqn = Integral(interior, u*u*u*(grad*v)*(grad*u), quad);
00098
00099 Expr bc = EssentialBC(northEdge, v*(u - pow(1.0 + sin(pi*x),0.25)),quad)
00100 + EssentialBC(southEdge + eastEdge + westEdge,
00101 v*(u - 1.0),quad);
00102
00103
00104 NonlinearProblem prob(mesh, eqn, bc, v, u, uNewt, vecType);
00105
00106 bool nonlinSolveOK = false;
00107 bool useNox = false;
00108 if (useNox)
00109 {
00110 nonlinSolveOK = noxNewtonSolve(prob, uNewt);
00111 }
00112 else
00113 {
00114 nonlinSolveOK = playaNewtonArmijoSolve(prob, uNewt);
00115 }
00116
00117 TEUCHOS_TEST_FOR_EXCEPT(!nonlinSolveOK);
00118
00119 FieldWriter writer = new VTKWriter("SteadyRadDiff2D");
00120 writer.addMesh(mesh);
00121 writer.addField("u", new ExprFieldWrapper(uNewt[0]));
00122 writer.write();
00123
00124 Expr uExact = pow(1.0 + sin(pi*x)*sinh(pi*y)/sinh(pi), 0.25);
00125 double err = L2Norm(mesh, interior, uNewt-uExact, quad);
00126 Out::root() << "error = " << setw(16) << err << endl;
00127
00128 double tol = 1.0e-4;
00129 Sundance::passFailTest(err, tol);
00130 }
00131 catch(std::exception& e)
00132 {
00133 Sundance::handleException(e);
00134 }
00135 Sundance::finalize(); return Sundance::testStatus();
00136 }
00137
00138
00139 bool inlineNewtonSolve(NonlinearProblem prob,
00140 Expr uNewt,
00141 int maxIters,
00142 double newtTol)
00143 {
00144 bool newtonConverged = false;
00145 LinearSolver<double> linSolver
00146 = LinearSolverBuilder::createSolver("amesos.xml");
00147
00148
00149 LinearOperator<double> J = prob.allocateJacobian();
00150 Vector<double> resid = J.range().createMember();
00151 Vector<double> newtonStep = J.domain().createMember();
00152
00153 for (int i=0; i<maxIters; i++)
00154 {
00155 prob.setEvalPoint(uNewt);
00156 prob.computeJacobianAndFunction(J, resid);
00157 SolverState<double> solveState
00158 = linSolver.solve(J, -1.0*resid, newtonStep);
00159
00160 TEUCHOS_TEST_FOR_EXCEPTION(solveState.finalState() != SolveConverged,
00161 std::runtime_error,
00162 "linear solve failed!");
00163
00164 addVecToDiscreteFunction(uNewt, newtonStep);
00165 double newtStepNorm = newtonStep.norm2();
00166 Out::root() << "|newt step| = " << newtStepNorm << endl;
00167 if (newtStepNorm < newtTol)
00168 {
00169 newtonConverged = true;
00170 break;
00171 }
00172 }
00173
00174 return newtonConverged;
00175 }
00176
00177 bool noxNewtonSolve(const NonlinearProblem& prob, Expr uNewt)
00178 {
00179
00180 ParameterXMLFileReader reader("nox-amesos.xml");
00181 ParameterList solverParams = reader.getParameters();
00182 NOXSolver solver(solverParams);
00183
00184 SolverState<double> stat = prob.solve(solver);
00185
00186 return stat.finalState()==SolveConverged;
00187 }
00188
00189 bool playaNewtonArmijoSolve(const NonlinearProblem& prob, Expr uNewt)
00190 {
00191
00192 LinearSolver<double> linSolver = LinearSolverBuilder::createSolver("amesos.xml");
00193 ParameterList solverParams("NewtonArmijoSolver");
00194 solverParams.set("Tau Relative", 1.0e-12);
00195 solverParams.set("Tau Absolute", 1.0e-12);
00196 solverParams.set("Alpha", 1.0e-4);
00197 solverParams.set("Verbosity", 3);
00198
00199 NonlinearSolver<double> solver = new NewtonArmijoSolver<double>(solverParams, linSolver);
00200
00201 SolverState<double> stat = prob.solve(solver);
00202 if (stat.finalState() != SolveConverged)
00203 {
00204 Out::os() << stat.finalMsg() << endl;
00205 }
00206
00207 return stat.finalState()==SolveConverged;
00208 }