00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "Sundance.hpp"
00032 #include "SundanceTransientStepProblem.hpp"
00033 #include "SundanceDoublingStepController.hpp"
00034 #include "SundanceProblemTesting.hpp"
00035
00036 Expr rhsF(const Expr& x, const Expr& t)
00037 {
00038 return -12*pow((1 - x)*cos(t) - x*cos(t),2)
00039 *pow(2 + (1 - x)*x*cos(t),2) +
00040 8*cos(t)*pow(2 + (1 - x)*x*cos(t),3) - (1 - x)*x*sin(t);
00041 }
00042
00043
00044 int main(int argc, char** argv)
00045 {
00046 try
00047 {
00048 int nx = 96;
00049
00050 Sundance::setOption("nx", nx, "Number of elements");
00051
00052 Sundance::init(&argc, &argv);
00053 int np = MPIComm::world().getNProc();
00054
00055
00056 VectorType<double> vecType = new EpetraVectorType();
00057
00058
00059
00060 MeshType meshType = new BasicSimplicialMeshType();
00061 MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx*np, meshType);
00062 Mesh mesh = mesher.getMesh();
00063
00064
00065
00066 CellFilter interior = new MaximalCellFilter();
00067 CellFilter points = new DimensionalCellFilter(0);
00068
00069 CellFilter right = points.coordSubset(0, 0.0);
00070 CellFilter left = points.coordSubset(0, 1.0);
00071 CellFilter midpoint = points.coordSubset(0, 0.5);
00072
00073
00074
00075 BasisFamily basis = new Lagrange(1);
00076 Expr u = new UnknownFunction(basis, "u");
00077 Expr v = new TestFunction(basis, "v");
00078
00079
00080
00081
00082 Expr dx = new Derivative(0);
00083 Expr x = new CoordExpr(0);
00084
00085
00086 QuadratureFamily quad = new GaussianQuadrature(4);
00087
00088
00089 Expr tPrev = new Sundance::Parameter(0.0);
00090 Expr t = new Sundance::Parameter(0.0);
00091 Expr dt = new Sundance::Parameter(0.0);
00092
00093 DiscreteSpace ds(mesh, basis, vecType);
00094 L2Projector proj(ds, 2.0+x*(1.0-x)*cos(t));
00095 Expr uPrev = proj.project();
00096 Expr uSoln = copyDiscreteFunction(uPrev);
00097
00098 Expr eqn = Integral(interior, v*(u-uPrev)
00099 + 2.0*dt*(dx*v)*(pow(u,3.0)*(dx*u) + pow(uPrev, 3.0)*(dx*uPrev))
00100 - 0.5*dt*v*(rhsF(x, tPrev) + rhsF(x, t)),
00101 quad);
00102 Expr bc = EssentialBC(left+right, v*(u+uPrev-4.0), quad);
00103
00104 NonlinearProblem stepProb(mesh, eqn, bc, v, u, uSoln, vecType);
00105
00106 TransientStepProblem prob(stepProb, tPrev, uPrev, t, uSoln, dt, 0);
00107
00108
00109 NonlinearSolver<double> solver
00110 = NonlinearSolverBuilder::createSolver("playa-newton-aztec-ml.xml");
00111
00112 double tol = 1.0e-4;
00113 double tFinal = 3.0*M_PI;
00114 StepControlParameters stepControl;
00115 stepControl.tStart_ = 0.0;
00116 stepControl.tStop_ = tFinal;
00117 stepControl.initialStepsize_ = 0.1;
00118 stepControl.maxStepsizeFactor_ = 100.0;
00119 stepControl.minStepsizeFactor_ = 1.0e-4;
00120 stepControl.verbosity_ = 4;
00121 stepControl.tau_ = tol;
00122
00123 FieldWriterFactory wf = new DSVWriterFactory();
00124 OutputControlParameters outputControl(wf, "ControlledTransient1D", tFinal/100.0, 0);
00125
00126 RCP<ExprComparisonBase> compare = rcp(new DefaultExprComparison());
00127
00128 DoublingStepController driver(prob, solver, stepControl, outputControl,
00129 compare);
00130
00131 driver.run();
00132
00133 Expr uFinal = copyDiscreteFunction(prob.uSoln());
00134 Expr uExact = 2.0 + x*(1.0-x)*cos(tFinal);
00135
00136 Expr errExpr = Integral(interior,
00137 pow(uFinal-uExact, 2),
00138 new GaussianQuadrature(2));
00139 double error = sqrt(evaluateIntegral(mesh, errExpr));
00140
00141 Out::root() << "error = " << error << endl;
00142
00143
00144 Sundance::passFailTest(error, tol);
00145 }
00146 catch(std::exception& e)
00147 {
00148 std::cerr << e.what() << std::endl;
00149 }
00150 Sundance::finalize(); return Sundance::testStatus();
00151 }
00152
00153
00154