libMesh
transient_ex1.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Transient Example 1 - Solving a Transient Linear System in Parallel</h1>
21 // \author Benjamin S. Kirk
22 // \date 2003
23 //
24 // This example shows how a simple, linear transient
25 // system can be solved in parallel. The system is simple
26 // scalar convection-diffusion with a specified external
27 // velocity. The initial condition is given, and the
28 // solution is advanced in time with a standard Crank-Nicolson
29 // time-stepping strategy.
30 
31 // C++ include files that we need
32 #include <iostream>
33 #include <algorithm>
34 #include <sstream>
35 #include <math.h>
36 
37 // Basic include file needed for the mesh functionality.
38 #include "libmesh/libmesh.h"
39 #include "libmesh/mesh.h"
40 #include "libmesh/mesh_refinement.h"
41 #include "libmesh/gmv_io.h"
42 #include "libmesh/equation_systems.h"
43 #include "libmesh/fe.h"
44 #include "libmesh/quadrature_gauss.h"
45 #include "libmesh/dof_map.h"
46 #include "libmesh/sparse_matrix.h"
47 #include "libmesh/numeric_vector.h"
48 #include "libmesh/dense_matrix.h"
49 #include "libmesh/dense_vector.h"
50 #include "libmesh/exodusII_io.h"
51 
52 // This example will solve a linear transient system,
53 // so we need to include the TransientLinearImplicitSystem definition.
54 #include "libmesh/linear_implicit_system.h"
55 #include "libmesh/transient_system.h"
56 #include "libmesh/vector_value.h"
57 
58 // The definition of a geometric element
59 #include "libmesh/elem.h"
60 
61 // Bring in everything from the libMesh namespace
62 using namespace libMesh;
63 
64 // Function prototype. This function will assemble the system
65 // matrix and right-hand-side at each time step. Note that
66 // since the system is linear we technically do not need to
67 // assemble the matrix at each time step, but we will anyway.
68 // In subsequent examples we will employ adaptive mesh refinement,
69 // and with a changing mesh it will be necessary to rebuild the
70 // system matrix.
71 void assemble_cd (EquationSystems & es,
72  const std::string & system_name);
73 
74 // Function prototype. This function will initialize the system.
75 // Initialization functions are optional for systems. They allow
76 // you to specify the initial values of the solution. If an
77 // initialization function is not provided then the default (0)
78 // solution is provided.
79 void init_cd (EquationSystems & es,
80  const std::string & system_name);
81 
82 // Exact solution function prototype. This gives the exact
83 // solution as a function of space and time. In this case the
84 // initial condition will be taken as the exact solution at time 0,
85 // as will the Dirichlet boundary conditions at time t.
86 Real exact_solution (const Real x,
87  const Real y,
88  const Real t);
89 
91  const Parameters & parameters,
92  const std::string &,
93  const std::string &)
94 {
95  return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
96 }
97 
98 
99 
100 // We can now begin the main program. Note that this
101 // example will fail if you are using complex numbers
102 // since it was designed to be run only with real numbers.
103 int main (int argc, char ** argv)
104 {
105  // Initialize libMesh.
106  LibMeshInit init (argc, argv);
107 
108  // This example requires a linear solver package.
109  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
110  "--enable-petsc, --enable-trilinos, or --enable-eigen");
111 
112  // This example requires Adaptive Mesh Refinement support - although
113  // it only refines uniformly, the refinement code used is the same
114  // underneath
115 #ifndef LIBMESH_ENABLE_AMR
116  libmesh_example_requires(false, "--enable-amr");
117 #else
118 
119  // Skip this 2D example if libMesh was compiled as 1D-only.
120  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
121 
122  // Read the mesh from file. This is the coarse mesh that will be used
123  // in example 10 to demonstrate adaptive mesh refinement. Here we will
124  // simply read it in and uniformly refine it 5 times before we compute
125  // with it.
126  //
127  // Create a mesh object, with dimension to be overridden later,
128  // distributed across the default MPI communicator.
129  Mesh mesh(init.comm());
130 
131  mesh.read ("mesh.xda");
132 
133  // Create a MeshRefinement object to handle refinement of our mesh.
134  // This class handles all the details of mesh refinement and coarsening.
135  MeshRefinement mesh_refinement (mesh);
136 
137  // Uniformly refine the mesh 5 times. This is the
138  // first time we use the mesh refinement capabilities
139  // of the library.
140  mesh_refinement.uniformly_refine (5);
141 
142  // Print information about the mesh to the screen.
143  mesh.print_info();
144 
145  // Create an equation systems object.
146  EquationSystems equation_systems (mesh);
147 
148  // Add a transient system to the EquationSystems
149  // object named "Convection-Diffusion".
151  equation_systems.add_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
152 
153  // Adds the variable "u" to "Convection-Diffusion". "u"
154  // will be approximated using first-order approximation.
155  system.add_variable ("u", FIRST);
156 
157  // Give the system a pointer to the matrix assembly
158  // and initialization functions.
159  system.attach_assemble_function (assemble_cd);
160  system.attach_init_function (init_cd);
161 
162  // Initialize the data structures for the equation system.
163  equation_systems.init ();
164 
165  // Prints information about the system to the screen.
166  equation_systems.print_info();
167 
168  // Write out the initial conditions.
169 #ifdef LIBMESH_HAVE_EXODUS_API
170  // If Exodus is available, we'll write all timesteps to the same file
171  // rather than one file per timestep.
172  std::string exodus_filename = "transient_ex1.e";
173  ExodusII_IO(mesh).write_equation_systems (exodus_filename, equation_systems);
174 #else
175  GMVIO(mesh).write_equation_systems ("out_000.gmv", equation_systems);
176 #endif
177 
178  // The Convection-Diffusion system requires that we specify
179  // the flow velocity. We will specify it as a RealVectorValue
180  // data type and then use the Parameters object to pass it to
181  // the assemble function.
182  equation_systems.parameters.set<RealVectorValue>("velocity") =
183  RealVectorValue (0.8, 0.8);
184 
185  // Solve the system "Convection-Diffusion". This will be done by
186  // looping over the specified time interval and calling the
187  // solve() member at each time step. This will assemble the
188  // system and call the linear solver.
189  const Real dt = 0.025;
190  system.time = 0.;
191 
192  for (unsigned int t_step = 0; t_step < 50; t_step++)
193  {
194  // Increment the time counter, set the time and the
195  // time step size as parameters in the EquationSystem.
196  system.time += dt;
197 
198  equation_systems.parameters.set<Real> ("time") = system.time;
199  equation_systems.parameters.set<Real> ("dt") = dt;
200 
201  // A pretty update message
202  libMesh::out << " Solving time step ";
203 
204  // Do fancy zero-padded formatting of the current time.
205  {
206  std::ostringstream out;
207 
208  out << std::setw(2)
209  << std::right
210  << t_step
211  << ", time="
212  << std::fixed
213  << std::setw(6)
214  << std::setprecision(3)
215  << std::setfill('0')
216  << std::left
217  << system.time
218  << "...";
219 
220  libMesh::out << out.str() << std::endl;
221  }
222 
223  // At this point we need to update the old
224  // solution vector. The old solution vector
225  // will be the current solution vector from the
226  // previous time step. We will do this by extracting the
227  // system from the EquationSystems object and using
228  // vector assignment. Since only TransientSystems
229  // (and systems derived from them) contain old solutions
230  // we need to specify the system type when we ask for it.
231  *system.old_local_solution = *system.current_local_solution;
232 
233  // Assemble & solve the linear system
234  equation_systems.get_system("Convection-Diffusion").solve();
235 
236  // Output every 10 timesteps to file.
237  if ((t_step+1)%10 == 0)
238  {
239 
240 #ifdef LIBMESH_HAVE_EXODUS_API
241  ExodusII_IO exo(mesh);
242  exo.append(true);
243  exo.write_timestep (exodus_filename, equation_systems, t_step+1, system.time);
244 #else
245  std::ostringstream file_name;
246 
247  file_name << "out_"
248  << std::setw(3)
249  << std::setfill('0')
250  << std::right
251  << t_step+1
252  << ".gmv";
253 
254 
255  GMVIO(mesh).write_equation_systems (file_name.str(),
256  equation_systems);
257 #endif
258  }
259  }
260 #endif // #ifdef LIBMESH_ENABLE_AMR
261 
262  // All done.
263  return 0;
264 }
265 
266 // We now define the function which provides the
267 // initialization routines for the "Convection-Diffusion"
268 // system. This handles things like setting initial
269 // conditions and boundary conditions.
271  const std::string & libmesh_dbg_var(system_name))
272 {
273  // It is a good idea to make sure we are initializing
274  // the proper system.
275  libmesh_assert_equal_to (system_name, "Convection-Diffusion");
276 
277  // Get a reference to the Convection-Diffusion system object.
279  es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
280 
281  // Project initial conditions at time 0
282  es.parameters.set<Real> ("time") = system.time = 0;
283 
284  system.project_solution(exact_value, libmesh_nullptr, es.parameters);
285 }
286 
287 
288 
289 // Now we define the assemble function which will be used
290 // by the EquationSystems object at each timestep to assemble
291 // the linear system for solution.
293  const std::string & system_name)
294 {
295  // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
296  libmesh_ignore(es);
297  libmesh_ignore(system_name);
298 
299 #ifdef LIBMESH_ENABLE_AMR
300  // It is a good idea to make sure we are assembling
301  // the proper system.
302  libmesh_assert_equal_to (system_name, "Convection-Diffusion");
303 
304  // Get a constant reference to the mesh object.
305  const MeshBase & mesh = es.get_mesh();
306 
307  // The dimension that we are running
308  const unsigned int dim = mesh.mesh_dimension();
309 
310  // Get a reference to the Convection-Diffusion system object.
312  es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
313 
314  // Get a constant reference to the Finite Element type
315  // for the first (and only) variable in the system.
316  FEType fe_type = system.variable_type(0);
317 
318  // Build a Finite Element object of the specified type. Since the
319  // FEBase::build() member dynamically creates memory we will
320  // store the object as a UniquePtr<FEBase>. This can be thought
321  // of as a pointer that will clean up after itself.
322  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
323  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
324 
325  // A Gauss quadrature rule for numerical integration.
326  // Let the FEType object decide what order rule is appropriate.
327  QGauss qrule (dim, fe_type.default_quadrature_order());
328  QGauss qface (dim-1, fe_type.default_quadrature_order());
329 
330  // Tell the finite element object to use our quadrature rule.
331  fe->attach_quadrature_rule (&qrule);
332  fe_face->attach_quadrature_rule (&qface);
333 
334  // Here we define some references to cell-specific data that
335  // will be used to assemble the linear system. We will start
336  // with the element Jacobian * quadrature weight at each integration point.
337  const std::vector<Real> & JxW = fe->get_JxW();
338  const std::vector<Real> & JxW_face = fe_face->get_JxW();
339 
340  // The element shape functions evaluated at the quadrature points.
341  const std::vector<std::vector<Real>> & phi = fe->get_phi();
342  const std::vector<std::vector<Real>> & psi = fe_face->get_phi();
343 
344  // The element shape function gradients evaluated at the quadrature
345  // points.
346  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
347 
348  // The XY locations of the quadrature points used for face integration
349  const std::vector<Point> & qface_points = fe_face->get_xyz();
350 
351  // A reference to the DofMap object for this system. The DofMap
352  // object handles the index translation from node and element numbers
353  // to degree of freedom numbers. We will talk more about the DofMap
354  // in future examples.
355  const DofMap & dof_map = system.get_dof_map();
356 
357  // Define data structures to contain the element matrix
358  // and right-hand-side vector contribution. Following
359  // basic finite element terminology we will denote these
360  // "Ke" and "Fe".
363 
364  // This vector will hold the degree of freedom indices for
365  // the element. These define where in the global system
366  // the element degrees of freedom get mapped.
367  std::vector<dof_id_type> dof_indices;
368 
369  // Here we extract the velocity & parameters that we put in the
370  // EquationSystems object.
371  const RealVectorValue velocity =
372  es.parameters.get<RealVectorValue> ("velocity");
373 
374  const Real dt = es.parameters.get<Real> ("dt");
375 
376  // Now we will loop over all the elements in the mesh that
377  // live on the local processor. We will compute the element
378  // matrix and right-hand-side contribution. Since the mesh
379  // will be refined we want to only consider the ACTIVE elements,
380  // hence we use a variant of the active_elem_iterator.
381  for (const auto & elem : mesh.active_local_element_ptr_range())
382  {
383  // Get the degree of freedom indices for the
384  // current element. These define where in the global
385  // matrix and right-hand-side this element will
386  // contribute to.
387  dof_map.dof_indices (elem, dof_indices);
388 
389  // Compute the element-specific data for the current
390  // element. This involves computing the location of the
391  // quadrature points (q_point) and the shape functions
392  // (phi, dphi) for the current element.
393  fe->reinit (elem);
394 
395  // Zero the element matrix and right-hand side before
396  // summing them. We use the resize member here because
397  // the number of degrees of freedom might have changed from
398  // the last element. Note that this will be the case if the
399  // element type is different (i.e. the last element was a
400  // triangle, now we are on a quadrilateral).
401  Ke.resize (dof_indices.size(),
402  dof_indices.size());
403 
404  Fe.resize (dof_indices.size());
405 
406  // Now we will build the element matrix and right-hand-side.
407  // Constructing the RHS requires the solution and its
408  // gradient from the previous timestep. This myst be
409  // calculated at each quadrature point by summing the
410  // solution degree-of-freedom values by the appropriate
411  // weight functions.
412  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
413  {
414  // Values to hold the old solution & its gradient.
415  Number u_old = 0.;
416  Gradient grad_u_old;
417 
418  // Compute the old solution & its gradient.
419  for (std::size_t l=0; l<phi.size(); l++)
420  {
421  u_old += phi[l][qp]*system.old_solution (dof_indices[l]);
422 
423  // This will work,
424  // grad_u_old += dphi[l][qp]*system.old_solution (dof_indices[l]);
425  // but we can do it without creating a temporary like this:
426  grad_u_old.add_scaled (dphi[l][qp], system.old_solution (dof_indices[l]));
427  }
428 
429  // Now compute the element matrix and RHS contributions.
430  for (std::size_t i=0; i<phi.size(); i++)
431  {
432  // The RHS contribution
433  Fe(i) += JxW[qp]*(
434  // Mass matrix term
435  u_old*phi[i][qp] +
436  -.5*dt*(
437  // Convection term
438  // (grad_u_old may be complex, so the
439  // order here is important!)
440  (grad_u_old*velocity)*phi[i][qp] +
441 
442  // Diffusion term
443  0.01*(grad_u_old*dphi[i][qp]))
444  );
445 
446  for (std::size_t j=0; j<phi.size(); j++)
447  {
448  // The matrix contribution
449  Ke(i,j) += JxW[qp]*(
450  // Mass-matrix
451  phi[i][qp]*phi[j][qp] +
452 
453  .5*dt*(
454  // Convection term
455  (velocity*dphi[j][qp])*phi[i][qp] +
456 
457  // Diffusion term
458  0.01*(dphi[i][qp]*dphi[j][qp]))
459  );
460  }
461  }
462  }
463 
464  // At this point the interior element integration has
465  // been completed. However, we have not yet addressed
466  // boundary conditions. For this example we will only
467  // consider simple Dirichlet boundary conditions imposed
468  // via the penalty method.
469  //
470  // The following loops over the sides of the element.
471  // If the element has no neighbor on a side then that
472  // side MUST live on a boundary of the domain.
473  {
474  // The penalty value.
475  const Real penalty = 1.e10;
476 
477  // The following loops over the sides of the element.
478  // If the element has no neighbor on a side then that
479  // side MUST live on a boundary of the domain.
480  for (auto s : elem->side_index_range())
481  if (elem->neighbor_ptr(s) == libmesh_nullptr)
482  {
483  fe_face->reinit(elem, s);
484 
485  for (unsigned int qp=0; qp<qface.n_points(); qp++)
486  {
487  const Number value = exact_solution (qface_points[qp](0),
488  qface_points[qp](1),
489  system.time);
490 
491  // RHS contribution
492  for (std::size_t i=0; i<psi.size(); i++)
493  Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
494 
495  // Matrix contribution
496  for (std::size_t i=0; i<psi.size(); i++)
497  for (std::size_t j=0; j<psi.size(); j++)
498  Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
499  }
500  }
501  }
502 
503  // If this assembly program were to be used on an adaptive mesh,
504  // we would have to apply any hanging node constraint equations
505  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
506 
507  // The element matrix and right-hand-side are now built
508  // for this element. Add them to the global matrix and
509  // right-hand-side vector. The SparseMatrix::add_matrix()
510  // and NumericVector::add_vector() members do this for us.
511  system.matrix->add_matrix (Ke, dof_indices);
512  system.rhs->add_vector (Fe, dof_indices);
513  }
514 
515  // That concludes the system matrix assembly routine.
516 #endif // #ifdef LIBMESH_ENABLE_AMR
517 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
UniquePtr< NumericVector< Number > > old_local_solution
All the values I need to compute my contribution to the simulation at hand.
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time)
Writes out the solution at a specific timestep.
Definition: exodusII_io.C:773
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
void add_scaled(const TypeVector< T2 > &, const T)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:624
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
int main(int argc, char **argv)
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
This class provides a specific system class.
Number exact_value(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
Definition: transient_ex1.C:90
const T & get(const std::string &) const
Definition: parameters.h:430
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
void assemble_cd(EquationSystems &es, const std::string &system_name)
void append(bool val)
If true, this flag will cause the ExodusII_IO object to attempt to open an existing file for writing...
Definition: exodusII_io.C:430
void init_cd(EquationSystems &es, const std::string &system_name)
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
PetscErrorCode Vec x
This is the MeshRefinement class.
std::string exodus_filename(unsigned number)
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
Real exact_solution(const Real x, const Real y, const Real t)
This is the exact solution that we are trying to obtain.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
void libmesh_ignore(const T &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
static const bool value
Definition: xdr_io.C:108
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
Number old_solution(const dof_id_type global_dof_number) const
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.