00001
00002
00003
00004 #include "Sundance.hpp"
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 Expr force(const double& epsilon, const Expr& x, const Expr& t)
00017 {
00018 const double pi = 4.0*atan(1.0);
00019
00020 Expr rtn = -8.0*epsilon*cos(2.0*pi*t)*
00021 pow(1.0 + pow(x,2.0)*epsilon*cos(2.0*pi*t),2.0)*
00022 (1.0 + 7.0*pow(x,2.0)*epsilon*cos(2.0*pi*t)) -
00023 2.0*pi*pow(x,2.0)*epsilon*sin(2.0*pi*t);
00024
00025 return rtn;
00026 }
00027
00028
00029 int main(int argc, char** argv)
00030 {
00031 const double pi = 4.0*atan(1.0);
00032
00033 try
00034 {
00035 Sundance::init(&argc, &argv);
00036 int np = MPIComm::world().getNProc();
00037
00038
00039 VectorType<double> vecType = new EpetraVectorType();
00040
00041
00042
00043 MeshType meshType = new BasicSimplicialMeshType();
00044 int nx = 32;
00045 MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx*np, meshType);
00046 Mesh mesh = mesher.getMesh();
00047
00048
00049
00050 CellFilter interior = new MaximalCellFilter();
00051 CellFilter points = new DimensionalCellFilter(0);
00052 CellFilter leftPoint = points.coordSubset(0, 0.0);
00053 CellFilter rightPoint = points.coordSubset(0, 1.0);
00054
00055
00056
00057
00058 BasisFamily bas = new Lagrange(1);
00059 Expr u = new UnknownFunction(bas, "u");
00060 Expr v = new TestFunction(bas, "v");
00061
00062
00063 Expr dx = new Derivative(0);
00064 Expr x = new CoordExpr(0);
00065
00066 double epsilon = 0.25;
00067
00068
00069
00070
00071
00072 Expr uStart = 1.0 + epsilon*x*x;
00073 DiscreteSpace discSpace(mesh, bas, vecType);
00074 L2Projector projector(discSpace, 1.0 + epsilon*x*x);
00075 Expr uPrev = projector.project();
00076
00077
00078 Expr uNewt = copyDiscreteFunction(uPrev, "uNewt");
00079
00080
00081 QuadratureFamily quad = new GaussianQuadrature(4);
00082
00083 int nSteps = nx;
00084 double dt = 1.0/((double) nSteps);
00085
00086
00087
00088 Expr t = new Sundance::Parameter(0.0);
00089 Expr tPrev = new Sundance::Parameter(0.0);
00090
00091
00092 Expr eqn = Integral(interior, v*(u-uPrev)
00093 + dt/2.0*(dx*v)*((dx*pow(u, 4.0))+(dx*pow(uPrev, 4.0)))
00094 - dt/2.0*v*(force(epsilon, x, t)+force(epsilon, x, tPrev)), quad);
00095
00096 Expr bc = EssentialBC(rightPoint, v*(u - 1.0 - epsilon*cos(2.0*pi*t)),quad);
00097
00098
00099 NonlinearProblem prob(mesh, eqn, bc, v, u, uNewt, vecType);
00100
00101 NonlinearSolver<double> solver
00102 = NonlinearSolverBuilder::createSolver("playa-newton-amesos.xml");
00103
00104
00105 FieldWriter w = new DSVWriter("transientNonlinear1D-0.dat");
00106 w.addMesh(mesh);
00107 w.addField("u", new ExprFieldWrapper(uPrev[0]));
00108 w.write();
00109
00110
00111 for (int i=0; i<nSteps; i++)
00112 {
00113
00114 Out::root() << "timestep #" << i << endl;
00115 t.setParameterValue((i+1)*dt);
00116 tPrev.setParameterValue(i*dt);
00117
00118 SolverState<double> state = prob.solve(solver);
00119
00120 TEUCHOS_TEST_FOR_EXCEPTION(state.finalState() != SolveConverged,
00121 std::runtime_error,
00122 "Nonlinear solve failed to converge: message=" << state.finalMsg());
00123
00124 updateDiscreteFunction(uNewt, uPrev);
00125
00126 FieldWriter writer = new MatlabWriter("transientNonlinear1D-"
00127 + Teuchos::toString(i+1) + ".dat");
00128 writer.addMesh(mesh);
00129 writer.addField("u", new ExprFieldWrapper(uPrev[0]));
00130 writer.write();
00131 }
00132
00133 double err = L2Norm(mesh, interior, uPrev - uStart, quad);
00134 Out::root() << "dt=" << setw(16) << dt << " error = " << setw(16) << err << endl;
00135
00136 double h = 1.0/(nx-1.0);
00137 double tol = 0.1*(dt*dt + h*h);
00138 Sundance::passFailTest(err, tol);
00139 }
00140 catch(std::exception& e)
00141 {
00142 Sundance::handleException(e);
00143 }
00144 Sundance::finalize(); return Sundance::testStatus();
00145 }