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 #include "Sundance.hpp"
00026 using Sundance::List;
00027
00028
00029 #if defined(HAVE_SUNDANCE_EXODUS) && defined(Trilinos_DATA_DIR)
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 int main(int argc, char** argv)
00040 {
00041 try
00042 {
00043 std::string meshFile="pipe2D-2";
00044 std::string outFile="TwoZoneTransport1";
00045 double K_f = 0.145;
00046 double c_f = 2.0;
00047 double rho_f = 0.865;
00048 double K_s = 0.61;
00049 double U0 = 20.0;
00050 double Q = 1.0;
00051 std::string heatSolverFile="aztec-ifpack.xml";
00052 std::string flowSolverFile="aztec-ml.xml";
00053
00054 Sundance::setOption("meshFile", meshFile, "mesh file (omit .exo suffix)");
00055 Sundance::setOption("outFile", outFile, "output file (omit .vtu suffix)");
00056 Sundance::setOption("K-fluid", K_f, "Thermal diffusivity of fluid");
00057 Sundance::setOption("rho-fluid", rho_f, "Density of fluid");
00058 Sundance::setOption("c-fluid", rho_f, "Specific heat of fluid");
00059 Sundance::setOption("K-solid", K_s, "Thermal diffusivity of solid");
00060 Sundance::setOption("U0", U0, "Average velocity of fluid");
00061 Sundance::setOption("Q", Q, "Heat generation rate in solid");
00062 Sundance::setOption("flow-solver", flowSolverFile, "flow solver XML file");
00063 Sundance::setOption("heat-solver", heatSolverFile, "heat solver XML file");
00064
00065 Sundance::init(&argc, &argv);
00066
00067 Out::root() << "K_f = " << K_f << endl;
00068 Out::root() << "rho_f = " << rho_f << endl;
00069 Out::root() << "c_f = " << c_f << endl;
00070 Out::root() << "K_s = " << K_s << endl;
00071 Out::root() << "Q = " << Q << endl;
00072 Out::root() << "Peclet number = " << rho_f*c_f*U0 / K_f << endl;
00073
00074
00075 VectorType<double> vecType = new EpetraVectorType();
00076
00077
00078 LinearSolver<double> flowSolver
00079 = LinearSolverBuilder::createSolver(flowSolverFile);
00080
00081 LinearSolver<double> heatSolver
00082 = LinearSolverBuilder::createSolver(heatSolverFile);
00083
00084
00085
00086 Out::root() << "reading mesh" << endl;
00087 MeshType meshType = new BasicSimplicialMeshType();
00088 MeshSource meshSrc
00089 = new ExodusMeshReader(meshFile, meshType);
00090 Mesh mesh = meshSrc.getMesh();
00091
00092
00093 CellFilter interior = new MaximalCellFilter();
00094 CellFilter fluid = interior.labeledSubset(1);
00095 CellFilter solid = interior.labeledSubset(2);
00096
00097 CellFilter sides = new DimensionalCellFilter(1);
00098
00099 CellFilter east = sides.labeledSubset(1);
00100 CellFilter west = sides.labeledSubset(2);
00101 CellFilter north = sides.labeledSubset(3);
00102 CellFilter south = sides.labeledSubset(4);
00103
00104
00105
00106 BasisFamily basis = new Lagrange(1);
00107 Expr phi = new UnknownFunction(basis, "phi");
00108 Expr v = new TestFunction(basis, "v");
00109
00110
00111 Expr grad = gradient(2);
00112
00113
00114 QuadratureFamily quad2 = new GaussianQuadrature(2);
00115 QuadratureFamily quad4 = new GaussianQuadrature(4);
00116
00117
00118 Expr flowEqn = Integral(fluid, (grad*phi)*(grad*v), quad2);
00119
00120 Expr h = new CellDiameterExpr();
00121 Expr x = new CoordExpr(0);
00122 Expr flowBC = EssentialBC(west, v*(phi-x*U0)/h, quad2) +
00123 EssentialBC(east, v*(phi-x*U0)/h, quad2);
00124
00125
00126 LinearProblem flowProb(mesh, flowEqn, flowBC, v, phi, vecType);
00127
00128 Out::root() << "solving flow" << endl;
00129 Expr phi0 = flowProb.solve(flowSolver);
00130
00131
00132 Expr T = new UnknownFunction(basis);
00133
00134 Expr n = CellNormalExpr(2, "n");
00135
00136 Expr heatEqn = Integral(solid, K_s*(grad*v)*(grad*T) - v*Q, quad2)
00137 + Integral(fluid,
00138 (grad*v)*(K_f*(grad*T) - c_f*rho_f*(grad*phi0)*T), quad2)
00139 + Integral(east, c_f*rho_f*v*T*n*(grad*phi0), quad2);
00140
00141 Expr heatBC = EssentialBC(west, v*T/h, quad2);
00142
00143
00144 LinearProblem heatProb(mesh, heatEqn, heatBC, v, T, vecType);
00145
00146 Out::root() << "solving heat" << endl;
00147 Expr T0 = heatProb.solve(heatSolver);
00148
00149 Out::root() << "writing output" << endl;
00150 FieldWriter w = new VTKWriter(outFile);
00151 w.addMesh(mesh);
00152 w.addField("phi", new ExprFieldWrapper(phi0[0]));
00153 w.addField("T", new ExprFieldWrapper(T0[0]));
00154 w.write();
00155
00156
00157
00158 Out::root() << "checking flux balance" << endl;
00159 Expr J = c_f*rho_f*T0*(grad*phi0) - K_f*(grad*T0);
00160 Expr fluxCheckNet = Integral(east+west+north+south, n*J, quad2);
00161 Expr sourceCheck = Integral(solid, Q, quad2);
00162
00163 double netFlux = evaluateIntegral(mesh, fluxCheckNet);
00164 double netSource = evaluateIntegral(mesh, sourceCheck);
00165
00166
00167 Out::root() << endl;
00168 Out::root() << "net heat flux out: " << netFlux << endl;
00169 Out::root() << "net heat production by source: " << netSource << endl;
00170 Out::root() << endl;
00171
00172
00173
00174
00175 Sundance::passFailTest(fabs(netFlux - netSource), 1.0e-2);
00176 }
00177 catch(std::exception& e)
00178 {
00179 Sundance::handleException(e);
00180 }
00181 Sundance::finalize(); return Sundance::testStatus();
00182
00183 return Sundance::testStatus();
00184 }
00185
00186 #else
00187
00188
00189 int main(int argc, char** argv)
00190 {
00191 Sundance::init(&argc, &argv);
00192 std::cout << "dummy PoissonDemo3D PASSED. Enable exodus to run the actual test" << std::endl;
00193 Sundance::finalize();
00194 return 0;
00195 }
00196
00197
00198 #endif
00199
00200