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, linearized by hand" << endl << endl;
00022
00023 VectorType<double> vecType = new EpetraVectorType();
00024
00025 MeshType meshType = new BasicSimplicialMeshType();
00026 MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx, meshType);
00027 Mesh mesh = mesher.getMesh();
00028
00029 CellFilter interior = new MaximalCellFilter();
00030 CellFilter sides = new DimensionalCellFilter(mesh.spatialDim()-1);
00031 CellFilter left = sides.subset(new CoordinateValueCellPredicate(0, 0.0));
00032 CellFilter right = sides.subset(new CoordinateValueCellPredicate(0, 1.0));
00033
00034 BasisFamily basis = new Lagrange(1);
00035 Expr w = new UnknownFunction(basis, "w");
00036 Expr v = new TestFunction(basis, "v");
00037
00038 Expr grad = gradient(1);
00039
00040 Expr x = new CoordExpr(0);
00041
00042
00043
00044 const double pi = 4.0*atan(1.0);
00045 Expr uExact = sin(pi*x);
00046 Expr R = pi*pi*uExact - lambda*exp(uExact);
00047
00048 QuadratureFamily quad4 = new GaussianQuadrature(4);
00049 QuadratureFamily quad2 = new GaussianQuadrature(2);
00050
00051 DiscreteSpace discSpace(mesh, basis, vecType);
00052 Expr uPrev = new DiscreteFunction(discSpace, 0.5);
00053 Expr stepVal = copyDiscreteFunction(uPrev);
00054
00055 Expr eqn
00056 = Integral(interior, (grad*v)*(grad*w) + (grad*v)*(grad*uPrev)
00057 - v*lambda*exp(uPrev)*(1.0+w) - v*R, quad4);
00058
00059 Expr h = new CellDiameterExpr();
00060 Expr bc = EssentialBC(left+right, v*(uPrev+w)/h, quad2);
00061
00062 LinearProblem prob(mesh, eqn, bc, v, w, vecType);
00063
00064 LinearSolver<double> linSolver
00065 = LinearSolverBuilder::createSolver("amesos.xml");
00066
00067 Out::root() << "Newton iteration" << endl;
00068 int maxIters = 20;
00069 Expr soln ;
00070 bool converged = false;
00071
00072 for (int i=0; i<maxIters; i++)
00073 {
00074
00075 prob.solve(linSolver, stepVal);
00076 Vector<double> stepVec = getDiscreteFunctionVector(stepVal);
00077 double deltaU = stepVec.norm2();
00078 Out::root() << "Iter=" << setw(3) << i << " ||Delta u||=" << setw(20)
00079 << deltaU << endl;
00080 addVecToDiscreteFunction(uPrev, stepVec);
00081 if (deltaU < convTol)
00082 {
00083 soln = uPrev;
00084 converged = true;
00085 break;
00086 }
00087 }
00088 TEUCHOS_TEST_FOR_EXCEPTION(!converged, std::runtime_error,
00089 "Newton iteration did not converge after "
00090 << maxIters << " iterations");
00091
00092 FieldWriter writer = new DSVWriter("HandCodedBratu.dat");
00093 writer.addMesh(mesh);
00094 writer.addField("soln", new ExprFieldWrapper(soln[0]));
00095 writer.write();
00096
00097 Out::root() << "Converged!" << endl << endl;
00098
00099 double L2Err = L2Norm(mesh, interior, soln-uExact, quad4);
00100 Out::root() << "L2 Norm of error: " << L2Err << endl;
00101
00102 Sundance::passFailTest(L2Err, 1.5/((double) nx*nx));
00103 }
00104 catch(exception& e)
00105 {
00106 Sundance::handleException(e);
00107 }
00108 Sundance::finalize();
00109 }
00110