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() << "Fixed-point iteration" << 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 u = new UnknownFunction(basis, "u");
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 uCur = copyDiscreteFunction(uPrev);
00054
00055 Expr eqn
00056 = Integral(interior, (grad*u)*(grad*v) - v*lambda*exp(uPrev) - v*R, quad4);
00057
00058 Expr h = new CellDiameterExpr();
00059 Expr bc = EssentialBC(left+right, v*u/h, quad4);
00060
00061 LinearProblem prob(mesh, eqn, bc, v, u, vecType);
00062
00063 Expr normSqExpr = Integral(interior, pow(u-uPrev, 2.0), quad2);
00064 Functional normSqFunc(mesh, normSqExpr, vecType);
00065 FunctionalEvaluator normSqEval = normSqFunc.evaluator(u, uCur);
00066
00067 LinearSolver<double> linSolver
00068 = LinearSolverBuilder::createSolver("amesos.xml");
00069
00070 Out::root() << "Fixed-point iteration" << endl;
00071 int maxIters = 20;
00072 Expr soln ;
00073 bool converged = false;
00074
00075 for (int i=0; i<maxIters; i++)
00076 {
00077
00078 prob.solve(linSolver, uCur);
00079
00080
00081 double deltaU = sqrt(normSqEval.evaluate());
00082 Out::root() << "Iter=" << setw(3) << i << " ||Delta u||=" << setw(20)
00083 << deltaU << endl;
00084
00085 if (deltaU < convTol)
00086 {
00087 soln = uCur;
00088 converged = true;
00089 break;
00090 }
00091
00092 Vector<double> uVec = getDiscreteFunctionVector(uCur);
00093
00094 setDiscreteFunctionVector(uPrev, uVec);
00095 }
00096 TEUCHOS_TEST_FOR_EXCEPTION(!converged, std::runtime_error,
00097 "Fixed point iteration did not converge after "
00098 << maxIters << " iterations");
00099
00100 FieldWriter writer = new DSVWriter("FixedPointBratu.dat");
00101 writer.addMesh(mesh);
00102 writer.addField("soln", new ExprFieldWrapper(soln[0]));
00103 writer.write();
00104
00105 Out::root() << "Converged!" << endl << endl;
00106
00107 double L2Err = L2Norm(mesh, interior, soln-uExact, quad4);
00108 Out::root() << "L2 Norm of error: " << L2Err << endl;
00109
00110 Sundance::passFailTest(L2Err, 1.5/((double) nx*nx));
00111 }
00112 catch(exception& e)
00113 {
00114 Sundance::handleException(e);
00115 }
00116 Sundance::finalize();
00117 }
00118