📄 ex18.php
字号:
<?php $root=""; ?><?php require($root."navigation.php"); ?><html><head> <?php load_style($root); ?></head> <body> <?php make_navigation("ex18",$root)?> <div class="content"><a name="comments"></a> <div class = "comment"><h1>Example 18 - Unsteady Navier-Stokes Equations with DiffSystem</h1><br><br>This example shows how the transient nonlinear problem fromexample 13 can be solved using the new (and experimental)DiffSystem class framework<br><br>Basic include files</div><div class ="fragment"><pre> #include "equation_systems.h" #include "error_vector.h" #include "getpot.h" #include "gmv_io.h" #include "kelly_error_estimator.h" #include "mesh.h" #include "mesh_generation.h" #include "mesh_refinement.h" #include "uniform_refinement_estimator.h" </pre></div><div class = "comment">Some (older) compilers do not offer full streamfunctionality, OStringStream works around this.</div><div class ="fragment"><pre> #include "o_string_stream.h" </pre></div><div class = "comment">The systems and solvers we may use</div><div class ="fragment"><pre> #include "naviersystem.h" #include "diff_solver.h" #include "euler_solver.h" #include "steady_solver.h" </pre></div><div class = "comment">The main program.</div><div class ="fragment"><pre> int main (int argc, char** argv) {</pre></div><div class = "comment">Initialize libMesh.</div><div class ="fragment"><pre> libMesh::init (argc, argv); #ifndef ENABLE_AMR std::cerr << "ERROR: This example requires libMesh to be\n" << "compiled with AMR support!" << std::endl; return 0; #else { </pre></div><div class = "comment">Parse the input file</div><div class ="fragment"><pre> GetPot infile("ex18.in"); </pre></div><div class = "comment">Read in parameters from the input file</div><div class ="fragment"><pre> const Real global_tolerance = infile("global_tolerance", 0.); const unsigned int nelem_target = infile("n_elements", 400); const bool transient = infile("transient", true); const Real deltat = infile("deltat", 0.005); unsigned int n_timesteps = infile("n_timesteps", 20); const unsigned int write_interval = infile("write_interval", 5); const unsigned int coarsegridsize = infile("coarsegridsize", 1); const unsigned int coarserefinements = infile("coarserefinements", 0); const unsigned int max_adaptivesteps = infile("max_adaptivesteps", 10); const unsigned int dim = infile("dimension", 2); assert (dim == 2 || dim == 3);</pre></div><div class = "comment">Create a n-dimensional mesh.</div><div class ="fragment"><pre> Mesh mesh (dim); </pre></div><div class = "comment">And an object to refine it</div><div class ="fragment"><pre> MeshRefinement mesh_refinement(mesh); mesh_refinement.coarsen_by_parents() = true; mesh_refinement.absolute_global_tolerance() = global_tolerance; mesh_refinement.nelem_target() = nelem_target; mesh_refinement.refine_fraction() = 0.3; mesh_refinement.coarsen_fraction() = 0.3; mesh_refinement.coarsen_threshold() = 0.1; </pre></div><div class = "comment">Use the MeshTools::Generation mesh generator to create a uniformgrid on the square [-1,1]^D. We instruct the mesh generatorto build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27elements in 3D. Building these higher-order elements allowsus to use higher-order approximation, as in example 3.</div><div class ="fragment"><pre> if (dim == 2) MeshTools::Generation::build_square (mesh, coarsegridsize, coarsegridsize, 0., 1., 0., 1., QUAD9); else if (dim == 3) MeshTools::Generation::build_cube (mesh, coarsegridsize, coarsegridsize, coarsegridsize, 0., 1., 0., 1., 0., 1., HEX27); mesh_refinement.uniformly_refine(coarserefinements); </pre></div><div class = "comment">Print information about the mesh to the screen.</div><div class ="fragment"><pre> mesh.print_info(); </pre></div><div class = "comment">Create an equation systems object.</div><div class ="fragment"><pre> EquationSystems equation_systems (mesh); </pre></div><div class = "comment">Declare the system "Navier-Stokes" and its variables.</div><div class ="fragment"><pre> NavierSystem & system = equation_systems.add_system<NavierSystem> ("Navier-Stokes"); </pre></div><div class = "comment">Solve this as a time-dependent or steady system</div><div class ="fragment"><pre> if (transient) system.time_solver = AutoPtr<TimeSolver>(new EulerSolver(system)); else { system.time_solver = AutoPtr<TimeSolver>(new SteadySolver(system)); assert(n_timesteps == 1); } </pre></div><div class = "comment">Initialize the system</div><div class ="fragment"><pre> equation_systems.init (); </pre></div><div class = "comment">Set the time stepping options</div><div class ="fragment"><pre> system.deltat = deltat; </pre></div><div class = "comment">And the nonlinear solver options</div><div class ="fragment"><pre> DiffSolver &solver = *(system.time_solver->diff_solver().get()); solver.quiet = infile("solver_quiet", true); solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15); solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3); solver.relative_residual_tolerance = infile("relative_residual_tolerance", 0.0); </pre></div><div class = "comment">And the linear solver options</div><div class ="fragment"><pre> solver.max_linear_iterations = infile("max_linear_iterations", 50000); solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3); </pre></div><div class = "comment">Print information about the system to the screen.</div><div class ="fragment"><pre> equation_systems.print_info(); </pre></div><div class = "comment">Now we begin the timestep loop to compute the time-accuratesolution of the equations.</div><div class ="fragment"><pre> for (unsigned int t_step=0; t_step != n_timesteps; ++t_step) {</pre></div><div class = "comment">A pretty update message</div><div class ="fragment"><pre> std::cout << " Solving time step " << t_step << ", time = " << system.time << std::endl; </pre></div><div class = "comment">Adaptively solve the timestep</div><div class ="fragment"><pre> unsigned int a_step = 0; for (; a_step != max_adaptivesteps; ++a_step) { system.solve(); ErrorVector error; AutoPtr<ErrorEstimator> error_estimator; </pre></div><div class = "comment">To solve to a tolerance in this problem weneed a better estimator than Kelly</div><div class ="fragment"><pre> if (global_tolerance != 0.) {</pre></div><div class = "comment">We can't adapt to both a tolerance and a meshsize at once</div><div class ="fragment"><pre> assert (nelem_target == 0); UniformRefinementEstimator *u = new UniformRefinementEstimator; </pre></div><div class = "comment">The lid-driven cavity problem isn't in H1, solets estimate H0 (i.e. L2) error</div><div class ="fragment"><pre> u->sobolev_order() = 0; error_estimator.reset(u); } else {</pre></div><div class = "comment">If we aren't adapting to a tolerance we need atarget mesh size</div><div class ="fragment"><pre> assert (nelem_target > 0); </pre></div><div class = "comment">Kelly is a lousy estimator to use for a problemnot in H1 - if we were doing more than a fewtimesteps we'd need to turn off or limit themaximum level of our adaptivity eventually</div><div class ="fragment"><pre> error_estimator.reset(new KellyErrorEstimator); } </pre></div><div class = "comment">Calculate error based on u and v but not p</div><div class ="fragment"><pre> error_estimator->component_scale.push_back(1.0); // u error_estimator->component_scale.push_back(1.0); // v if (dim == 3) error_estimator->component_scale.push_back(1.0); // w error_estimator->component_scale.push_back(0.0); // p error_estimator->estimate_error(system, error); ⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -