00001 #include "Sundance.hpp"
00002
00003 #if defined(HAVE_SUNDANCE_EXODUS) && defined(Trilinos_DATA_DIR)
00004
00005
00006 Array<double> runit(const string& meshName)
00007 {
00008 VectorType<double> vecType = new EpetraVectorType();
00009 MeshType meshType = new BasicSimplicialMeshType();
00010 MeshSource mesher = new ExodusMeshReader(meshName, meshType);
00011 Mesh mesh = mesher.getMesh();
00012
00013 BasisFamily P2 = new Lagrange(2);
00014 BasisFamily P1 = new Lagrange(1);
00015 QuadratureFamily quad = new GaussianQuadrature(4);
00016
00017 Expr ux = new UnknownFunction(P2);
00018 Expr uy = new UnknownFunction(P2);
00019 Expr vx = new TestFunction(P2);
00020 Expr vy = new TestFunction(P2);
00021 Expr u = List(ux, uy);
00022 Expr v = List(vx, vy);
00023
00024 Expr p = new UnknownFunction(P1);
00025 Expr q = new TestFunction(P1);
00026
00027 CellFilter interior = new MaximalCellFilter();
00028 CellFilter sides = new DimensionalCellFilter(1);
00029 CellFilter pts = new DimensionalCellFilter(0);
00030 CellFilter peg = pts.labeledSubset(1);
00031 CellFilter inner = sides.labeledSubset(2);
00032 CellFilter outer = sides.labeledSubset(1);
00033
00034 Expr grad = gradient(2);
00035 Expr x = new CoordExpr(0);
00036 Expr y = new CoordExpr(1);
00037 Expr r = sqrt(x*x+y*y);
00038 double uIn = 0.25;
00039 double uOut = 1.0;
00040 Expr uInnerX = -uIn*y/r;
00041 Expr uInnerY = uIn*x/r;
00042 Expr uOuterX = -uOut*y/r;
00043 Expr uOuterY = uOut*x/r;
00044 Expr uInner = List(uInnerX, uInnerY);
00045 Expr uOuter = List(uOuterX, uOuterY);
00046
00047 Expr nu = 1.0;
00048 double C1 = 20.0;
00049 double C2 = 20.0;
00050
00051
00052 Expr eqn = Integral(interior,
00053 nu * colonProduct( outerProduct(grad, u), outerProduct(grad, v) )
00054 - p*div(v) - q*div(u), quad)
00055 + NitscheStokesNoSlipBC(inner, quad, nu, v, q, u, p, uInner, C1, C2)
00056 + NitscheStokesNoSlipBC(outer, quad, nu, v, q, u, p, uOuter, C1, C2);
00057
00058 Expr dummyBC;
00059
00060
00061 Expr MpEqn = Integral(interior, p*q, quad);
00062 Expr ApEqn = Integral(interior, (grad*p)*(grad*q), quad)
00063 + Integral(peg, p*q, quad);
00064 Expr FpEqn = Integral(interior, (grad*p)*(grad*q), quad)
00065 + Integral(peg, p*q, quad);
00066
00067
00068 Array<Block> unkBlocks = tuple(Block(u, vecType), Block(p, vecType));
00069 Array<Block> varBlocks = tuple(Block(v, vecType), Block(q, vecType));
00070
00071
00072 LinearProblem prob( mesh , eqn , dummyBC , varBlocks, unkBlocks );
00073
00074
00075 LinearProblem MpProb(mesh, MpEqn, dummyBC, q, p, vecType);
00076 LinearProblem ApProb(mesh, ApEqn, dummyBC, q, p, vecType);
00077 LinearProblem FpProb(mesh, FpEqn, dummyBC, q, p, vecType);
00078
00079
00080
00081 ParameterXMLFileReader outerSolverReader("pcd-outer-belos.xml");
00082 ParameterList outerSolverParams = outerSolverReader.getParameters();
00083 ParameterList precParams = outerSolverParams.sublist("PCD Preconditioner");
00084 LinearSolver<double> outerSolver = LinearSolverBuilder::createSolver(outerSolverParams);
00085
00086
00087 PreconditionerFactory<double> outerPrec
00088 = new PCDPreconditionerFactory(precParams, MpProb, ApProb, FpProb);
00089 outerSolver.setUserPrec(outerPrec);
00090
00091
00092 Expr soln = prob.solve(outerSolver);
00093
00094
00095
00096 double c1 = -2.0/3.0*(uIn-2.0*uOut);
00097 double c2 = 1.0/3.0*(2.0*uIn - uOut);
00098 Expr uMagExact = c1*r + c2/r;
00099 Expr uxExact = -uMagExact*y/r;
00100 Expr uyExact = uMagExact*x/r;
00101 Expr uExact = List(uxExact, uyExact);
00102
00103 double area = L2Norm(mesh, interior, 1.0, quad);
00104 double uErrorNorm = L2Norm(mesh, interior, soln[0]-uExact, quad)/area;
00105 double pErrorNorm = L2Norm(mesh, interior, soln[1], quad)/area;
00106
00107
00108 Expr h = new CellDiameterExpr();
00109 double hAvg = L2Norm(mesh, interior, h, quad)/area;
00110
00111
00112 Out::root() << "average h = " << hAvg << std::endl;
00113 Out::root() << "velocity error norm = " << uErrorNorm << std::endl;
00114 Out::root() << "pressure error norm = " << pErrorNorm << std::endl;
00115
00116
00117
00118
00119
00120 DiscreteSpace P1Space( mesh , P1 , vecType );
00121 L2Projector proj_ux( P1Space , soln[0][0] );
00122 L2Projector proj_uy( P1Space , soln[0][1] );
00123 L2Projector proj_uxEx( P1Space , uxExact );
00124 L2Projector proj_uyEx( P1Space , uyExact );
00125
00126 Expr ux_P1 = proj_ux.project( );
00127 Expr uy_P1 = proj_uy.project( );
00128
00129 Expr uxEx_P1 = proj_uxEx.project( );
00130 Expr uyEx_P1 = proj_uyEx.project( );
00131
00132 FieldWriter w = new VTKWriter("PCDStokes-" + meshName);
00133 w.addMesh(mesh);
00134 w.addField( "ux" , new ExprFieldWrapper( ux_P1 ) );
00135 w.addField( "uy" , new ExprFieldWrapper( uy_P1 ) );
00136 w.addField( "p" , new ExprFieldWrapper( soln[1] ) );
00137 w.addField( "uxEx" , new ExprFieldWrapper( uxEx_P1 ) );
00138 w.addField( "uyEx" , new ExprFieldWrapper( uyEx_P1 ) );
00139 w.write();
00140
00141 return tuple<double>(hAvg, uErrorNorm, pErrorNorm);
00142 }
00143
00144
00145
00146
00147
00148
00149
00150
00151 int main( int argc , char **argv )
00152 {
00153 try
00154 {
00155 Sundance::init( &argc , &argv );
00156
00157 Array<Array<double> > results;
00158 for (int n=0; n<=1; n++)
00159 {
00160 string meshName = "diskWithHole2D-" + Teuchos::toString(n);
00161 Array<double> x = runit(meshName);
00162 results.append(x);
00163 }
00164
00165 for (int n=0; n<results.size(); n++)
00166 {
00167 Out::root() << setw(12) << results[n][0] << setw(20) << results[n][1]
00168 << setw(20) << results[n][2] << endl;
00169 }
00170
00171 Sundance::passFailTest(results[0][1], 0.001);
00172 }
00173 catch (exception &e) {
00174 Out::os() << e.what() << std::endl;
00175 }
00176
00177 Sundance::finalize();
00178 return Sundance::testStatus();
00179 }
00180
00181
00182
00183
00184 #else
00185
00186
00187 int main(int argc, char** argv)
00188 {
00189 Sundance::init(&argc, &argv);
00190 std::cout << "dummy PCDStokes PASSED. Enable exodus to run the actual test" << std::endl;
00191 Sundance::finalize();
00192 return 0;
00193 }
00194
00195
00196 #endif