00001 #include "Sundance.hpp"
00002
00003
00004
00005
00006
00007 int main(int argc, char** argv)
00008 {
00009 try
00010 {
00011 int nx = 32;
00012 double convTol = 1.0e-8;
00013 double lambda = 0.5;
00014 Sundance::setOption("nx", nx, "Number of elements");
00015 Sundance::setOption("tol", convTol, "Convergence tolerance");
00016 Sundance::setOption("lambda", lambda, "Lambda (parameter in Bratu's equation)");
00017
00018 Sundance::init(&argc, &argv);
00019
00020 Out::root() << "Bratu problem (lambda=" << lambda << ")" << endl;
00021 Out::root() << "Newton's method with automated linearization"
00022 << endl << endl;
00023
00024 VectorType<double> vecType = new EpetraVectorType();
00025
00026 MeshType meshType = new BasicSimplicialMeshType();
00027 MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx, meshType);
00028 Mesh mesh = mesher.getMesh();
00029
00030 CellFilter interior = new MaximalCellFilter();
00031 CellFilter sides = new DimensionalCellFilter(mesh.spatialDim()-1);
00032 CellFilter left = sides.subset(new CoordinateValueCellPredicate(0, 0.0));
00033 CellFilter right = sides.subset(new CoordinateValueCellPredicate(0, 1.0));
00034
00035 BasisFamily basis = new Lagrange(1);
00036 Expr u = new UnknownFunction(basis, "w");
00037 Expr v = new TestFunction(basis, "v");
00038
00039 Expr grad = gradient(1);
00040
00041 Expr x = new CoordExpr(0);
00042
00043 const double pi = 4.0*atan(1.0);
00044 Expr uExact = sin(pi*x);
00045 Expr R = pi*pi*uExact - lambda*exp(uExact);
00046
00047 QuadratureFamily quad4 = new GaussianQuadrature(4);
00048 QuadratureFamily quad2 = new GaussianQuadrature(2);
00049
00050 DiscreteSpace discSpace(mesh, basis, vecType);
00051 Expr uPrev = new DiscreteFunction(discSpace, 0.5);
00052
00053 Expr eqn
00054 = Integral(interior, (grad*v)*(grad*u) - v*lambda*exp(u) - v*R, quad4);
00055
00056 Expr h = new CellDiameterExpr();
00057 Expr bc = EssentialBC(left+right, v*u/h, quad2);
00058
00059 NonlinearProblem prob(mesh, eqn, bc, v, u, uPrev, vecType);
00060
00061 LinearSolver<double> linSolver
00062 = LinearSolverBuilder::createSolver("amesos.xml");
00063
00064 Out::root() << "Newton iteration" << endl;
00065 int maxIters = 20;
00066 Expr soln ;
00067 bool converged = false;
00068
00069 LinearOperator<double> J = prob.allocateJacobian();
00070 Vector<double> residVec = J.range().createMember();
00071 Vector<double> stepVec;
00072
00073 for (int i=0; i<maxIters; i++)
00074 {
00075 prob.setInitialGuess(uPrev);
00076 prob.computeJacobianAndFunction(J, residVec);
00077
00078 linSolver.solve(J, -1.0*residVec, stepVec);
00079
00080 double deltaU = stepVec.norm2();
00081 Out::root() << "Iter=" << setw(3) << i << " ||Delta u||=" << setw(20)
00082 << deltaU << endl;
00083 addVecToDiscreteFunction(uPrev, stepVec);
00084 if (deltaU < convTol)
00085 {
00086 soln = uPrev;
00087 converged = true;
00088 break;
00089 }
00090 }
00091 TEUCHOS_TEST_FOR_EXCEPTION(!converged, std::runtime_error,
00092 "Newton iteration did not converge after "
00093 << maxIters << " iterations");
00094
00095 FieldWriter writer = new DSVWriter("AutoLinearizedBratu.dat");
00096 writer.addMesh(mesh);
00097 writer.addField("soln", new ExprFieldWrapper(soln[0]));
00098 writer.write();
00099
00100 Out::root() << "Converged!" << endl << endl;
00101
00102 double L2Err = L2Norm(mesh, interior, soln-uExact, quad4);
00103 Out::root() << "L2 Norm of error: " << L2Err << endl;
00104
00105 Sundance::passFailTest(L2Err, 1.5/((double) nx*nx));
00106 }
00107 catch(exception& e)
00108 {
00109 Sundance::handleException(e);
00110 }
00111 Sundance::finalize();
00112 }
00113
