📄 ex10.php
字号:
<?php $root=""; ?><?php require($root."navigation.php"); ?><html><head> <?php load_style($root); ?></head> <body> <?php make_navigation("ex10",$root)?> <div class="content"><a name="comments"></a> <div class = "comment"><h1>Example 10 - Solving a Transient System with Adaptive Mesh Refinement</h1><br><br>This example shows how a simple, linear transientsystem can be solved in parallel. The system is simplescalar convection-diffusion with a specified externalvelocity. The initial condition is given, and thesolution is advanced in time with a standard Crank-Nicholsontime-stepping strategy.<br><br>Also, we use this example to demonstrate writing out and readingin of solutions. We do 25 time steps, then save the solutionand do another 25 time steps starting from the saved solution. <br><br>C++ include files that we need</div><div class ="fragment"><pre> #include <iostream> #include <algorithm> #include <cmath> </pre></div><div class = "comment">Basic include file needed for the mesh functionality.</div><div class ="fragment"><pre> #include "libmesh.h" #include "mesh.h" #include "mesh_refinement.h" #include "gmv_io.h" #include "equation_systems.h" #include "fe.h" #include "quadrature_gauss.h" #include "dof_map.h" #include "sparse_matrix.h" #include "numeric_vector.h" #include "dense_matrix.h" #include "dense_vector.h" #include "getpot.h" </pre></div><div class = "comment">Some (older) compilers do not offer full stream functionality, \p OStringStream works around this.Check example 9 for details.</div><div class ="fragment"><pre> #include "o_string_stream.h" </pre></div><div class = "comment">This example will solve a linear transient system,so we need to include the \p TransientLinearImplicitSystem definition.</div><div class ="fragment"><pre> #include "transient_system.h" #include "linear_implicit_system.h" #include "vector_value.h" </pre></div><div class = "comment">To refine the mesh we need an \p ErrorEstimatorobject to figure out which elements to refine.</div><div class ="fragment"><pre> #include "error_vector.h" #include "kelly_error_estimator.h" </pre></div><div class = "comment">The definition of a geometric element</div><div class ="fragment"><pre> #include "elem.h" </pre></div><div class = "comment">Function prototype. This function will assemble the systemmatrix and right-hand-side at each time step. Note thatsince the system is linear we technically do not need toassmeble the matrix at each time step, but we will anyway.In subsequent examples we will employ adaptive mesh refinement,and with a changing mesh it will be necessary to rebuild thesystem matrix.</div><div class ="fragment"><pre> void assemble_cd (EquationSystems& es, const std::string& system_name); </pre></div><div class = "comment">Function prototype. This function will initialize the system.Initialization functions are optional for systems. They allowyou to specify the initial values of the solution. If aninitialization function is not provided then the default (0)solution is provided.</div><div class ="fragment"><pre> void init_cd (EquationSystems& es, const std::string& system_name); </pre></div><div class = "comment">Exact solution function prototype. This gives the exactsolution as a function of space and time. In this case theinitial condition will be taken as the exact solution at time 0,as will the Dirichlet boundary conditions at time t.</div><div class ="fragment"><pre> Real exact_solution (const Real x, const Real y, const Real t); Number exact_value (const Point& p, const Parameters& parameters, const std::string&, const std::string&) { return exact_solution(p(0), p(1), parameters.get<Real> ("time")); } </pre></div><div class = "comment">Begin the main program. Note that the firststatement in the program throws an error ifyou are in complex number mode, since thisexample is only intended to work with realnumbers.</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">Brief message to the user regarding the program nameand command line arguments.<br><br>Use commandline parameter to specify if we are toread in an initial solution or generate it ourself</div><div class ="fragment"><pre> std::cout << "Usage:\n" <<"\t " << argv[0] << " -init_timestep 0\n" << "OR\n" <<"\t " << argv[0] << " -read_solution -init_timestep 26\n" << std::endl; std::cout << "Running: " << argv[0]; for (int i=1; i<argc; i++) std::cout << " " << argv[i]; std::cout << std::endl << std::endl; </pre></div><div class = "comment">Create a GetPot object to parse the command line</div><div class ="fragment"><pre> GetPot command_line (argc, argv); </pre></div><div class = "comment">This boolean value is obtained from the command line, it is trueif the flag "-read_solution" is present, false otherwise.It indicates whether we are going to read inthe mesh and solution files "saved_mesh.xda" and "saved_solution.xda"or whether we are going to start from scratch by just reading"mesh.xda"</div><div class ="fragment"><pre> const bool read_solution = command_line.search("-read_solution"); </pre></div><div class = "comment">This value is also obtained from the commandline and it specifies theinitial value for the t_step looping variable. We mustdistinguish between the two cases here, whether we read in the solution or we started from scratch, so that we do not overwrite thegmv output files.</div><div class ="fragment"><pre> unsigned int init_timestep = 0; </pre></div><div class = "comment">Search the command line for the "init_timestep" flag and if it ispresent, set init_timestep accordingly.</div><div class ="fragment"><pre> if(command_line.search("-init_timestep")) init_timestep = command_line.next(0); else { std::cerr << "ERROR: Initial timestep not specified\n" << std::endl; </pre></div><div class = "comment">This handy function will print the file name, line number,and then abort. Currrently the library does not use C++exception handling.</div><div class ="fragment"><pre> error(); } </pre></div><div class = "comment">This value is also obtained from the command line, and specifiesthe number of time steps to take.</div><div class ="fragment"><pre> unsigned int n_timesteps = 0; </pre></div><div class = "comment">Again do a search on the command line for the argument</div><div class ="fragment"><pre> if(command_line.search("-n_timesteps")) n_timesteps = command_line.next(0); else { std::cerr << "ERROR: Number of timesteps not specified\n" << std::endl; error(); } </pre></div><div class = "comment">Create a two-dimensional mesh.</div><div class ="fragment"><pre> Mesh mesh (2); </pre></div><div class = "comment">Create an equation systems object.</div><div class ="fragment"><pre> EquationSystems equation_systems (mesh); MeshRefinement mesh_refinement (mesh); </pre></div><div class = "comment">First we process the case where we do not read in the solution</div><div class ="fragment"><pre> if(!read_solution) {</pre></div><div class = "comment">Read the mesh from file.</div><div class ="fragment"><pre> mesh.read ("mesh.xda"); </pre></div><div class = "comment">Uniformly refine the mesh 5 times</div><div class ="fragment"><pre> if(!read_solution) mesh_refinement.uniformly_refine (5); </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">Declare the system and its variables.Begin by creating a transient systemnamed "Convection-Diffusion".</div><div class ="fragment"><pre> TransientLinearImplicitSystem & system = equation_systems.add_system<TransientLinearImplicitSystem> ("Convection-Diffusion"); </pre></div><div class = "comment">Adds the variable "u" to "Convection-Diffusion". "u"will be approximated using first-order approximation.</div><div class ="fragment"><pre> system.add_variable ("u", FIRST); </pre></div><div class = "comment">Give the system a pointer to the matrix assemblyand initialization functions.</div><div class ="fragment"><pre> system.attach_assemble_function (assemble_cd); system.attach_init_function (init_cd); </pre></div><div class = "comment">Initialize the data structures for the equation system.</div><div class ="fragment"><pre> equation_systems.init (); } </pre></div><div class = "comment">Otherwise we read in the solution and mesh</div><div class ="fragment"><pre> else {</pre></div><div class = "comment">Read in the mesh stored in "saved_mesh.xda"</div><div class ="fragment"><pre> mesh.read("saved_mesh.xda"); </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">Read in the solution stored in "saved_solution.xda"</div><div class ="fragment"><pre> equation_systems.read("saved_solution.xda", libMeshEnums::READ); </pre></div><div class = "comment">Get a reference to the system so that we can call update() on it</div><div class ="fragment"><pre> TransientLinearImplicitSystem & system = equation_systems.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion"); </pre></div><div class = "comment">We need to call update to put system in a consistent statewith the solution that was read in</div><div class ="fragment"><pre> system.update(); </pre></div><div class = "comment">Attach the same matrix assembly function as above. Note, we do nothave to attach an init() function since we are initializing the
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -