00001 #include "Sundance.hpp"
00002 #include "PlayaNonlinearSolverBuilder.hpp"
00003
00004
00005
00006
00007
00008 int main(int argc, char** argv)
00009 {
00010 try
00011 {
00012 int nx = 32;
00013 double convTol = 1.0e-8;
00014 int nSteps = 5;
00015 double lambdaMax = 0.5;
00016 Sundance::setOption("nx", nx, "Number of elements");
00017 Sundance::setOption("tol", convTol, "Convergence tolerance");
00018 Sundance::setOption("lambda-max", lambdaMax,
00019 "final lambda (parameter in Bratu's equation)");
00020 Sundance::setOption("nSteps", nSteps,
00021 "number of steps in lambda (continuation from lambda=0 to lambda=lambdaMax)");
00022
00023 Sundance::init(&argc, &argv);
00024
00025 Out::root() << "Bratu problem with continuation (lambda=[0, " << lambdaMax << "] in "
00026 << nSteps << " steps)" << endl;
00027 Out::root() << "Newton's method with automated linearization"
00028 << endl << endl;
00029
00030 VectorType<double> vecType = new EpetraVectorType();
00031
00032 MeshType meshType = new BasicSimplicialMeshType();
00033 MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx, meshType);
00034 Mesh mesh = mesher.getMesh();
00035
00036 CellFilter interior = new MaximalCellFilter();
00037 CellFilter sides = new DimensionalCellFilter(mesh.spatialDim()-1);
00038 CellFilter left = sides.subset(new CoordinateValueCellPredicate(0, 0.0));
00039 CellFilter right = sides.subset(new CoordinateValueCellPredicate(0, 1.0));
00040
00041 BasisFamily basis = new Lagrange(1);
00042 Expr u = new UnknownFunction(basis, "w");
00043 Expr v = new TestFunction(basis, "v");
00044
00045 Expr grad = gradient(1);
00046
00047 Expr x = new CoordExpr(0);
00048
00049 Expr lambda = new Sundance::Parameter(0.0);
00050 const double pi = 4.0*atan(1.0);
00051 Expr uExact = sin(pi*x);
00052 Expr R = pi*pi*uExact - lambda*exp(uExact);
00053
00054 QuadratureFamily quad4 = new GaussianQuadrature(4);
00055 QuadratureFamily quad2 = new GaussianQuadrature(2);
00056
00057 DiscreteSpace discSpace(mesh, basis, vecType);
00058 Expr uPrev = new DiscreteFunction(discSpace, 0.5);
00059
00060 Expr eqn
00061 = Integral(interior, (grad*v)*(grad*u) - v*lambda*exp(u) - v*R, quad4);
00062
00063 Expr h = new CellDiameterExpr();
00064 Expr bc = EssentialBC(left+right, v*u/h, quad2);
00065
00066 NonlinearProblem prob(mesh, eqn, bc, v, u, uPrev, vecType);
00067
00068 NonlinearSolver<double> solver
00069 = NonlinearSolverBuilder::createSolver("playa-newton-amesos.xml");
00070
00071 Expr soln = uPrev;
00072 for (int n=0; n<nSteps; n++)
00073 {
00074 double lambdaVal = n*lambdaMax/(nSteps-1.0);
00075
00076 lambda.setParameterValue(lambdaVal);
00077 Out::root() << "continuation step n=" << n
00078 << " of " << nSteps << ", lambda="
00079 << lambdaVal << endl;
00080
00081 SolverState<double> state = prob.solve(solver);
00082
00083 TEUCHOS_TEST_FOR_EXCEPTION(state.finalState() != SolveConverged,
00084 std::runtime_error,
00085 "Nonlinear solve failed to converge: message=" << state.finalMsg());
00086
00087 Expr soln = uPrev;
00088 FieldWriter writer = new DSVWriter("ContinuationBratu-"
00089 + Teuchos::toString(n) + ".dat");
00090 writer.addMesh(mesh);
00091 writer.addField("soln", new ExprFieldWrapper(soln[0]));
00092 writer.write();
00093 }
00094
00095
00096 double L2Err = L2Norm(mesh, interior, soln-uExact, quad4);
00097 Out::root() << "L2 Norm of error: " << L2Err << endl;
00098
00099 Sundance::passFailTest(L2Err, 1.5/((double) nx*nx));
00100 }
00101 catch(std::exception& e)
00102 {
00103 Sundance::handleException(e);
00104 }
00105 Sundance::finalize();
00106 }
00107