libMesh
adjoints_ex6.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 
15 
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20 // <h1>Adjoints Example 6 - Poisson Equation in square domain with AdjointRefinementErrorEstimator</h1>
21 // \author Vikram Garg
22 // \date 2017
23 //
24 // This example solves the Poisson equation, whose solution displays a
25 // sharp layer, with QoI based adjoint error estimation and adaptive
26 // mesh refinement. The exact QoI value is in poisson.in. This example
27 // also illustrates the use of the adjoint Dirichlet boundary
28 // condition capability, necessary for handling flux QoIs. We access
29 // the adjoint capabilities of libMesh via the DiffSystem
30 // framework. This file (adjoints_ex6.C) contains the declaration of
31 // mesh and equation system objects, poissonsystem.C contains the
32 // assembly of the system. Postprocessing to compute the QoI is done
33 // in element_postprocess.C. There is no need for element and side
34 // qoi derivative functions, since the adjoint RHS is supplied by the
35 // adjoint dirichlet boundary condition.
36 
37 // WARNING: Adjoint-dirichlet based weighted flux quantities of
38 // interest are computed internally by FEMSystem as R^h(u^h, L^h)
39 // where R^h is the physics residual, u^h is the solution, and L^h the
40 // lift function at the current mesh. This is appropriate for
41 // unstabilized weak residuals. Non-FEMSystem users and stabilized
42 // method users will need to override the default flux QoI behavior.
43 
44 // An input file named "general.in" is provided which allows the user
45 // to set several parameters for the solution so that the problem can
46 // be re-run without a re-compile. The solution technique employed is
47 // to have a refinement loop with a linear (forward and adjoint) solve
48 // inside followed by a refinement of the grid and projection of the
49 // solution to the new grid In the final loop iteration, there is no
50 // additional refinement after the solve. In the input file
51 // "general.in", the variable "max_adaptivesteps" controls the number
52 // of refinement steps, and "refine_fraction" / "coarsen_fraction"
53 // determine the number of elements which will be refined / coarsened
54 // at each step.
55 
56 // C++ includes
57 #include <iostream>
58 #include <iomanip>
59 
60 // General libMesh includes
61 #include "libmesh/equation_systems.h"
62 #include "libmesh/linear_solver.h"
63 #include "libmesh/error_vector.h"
64 #include "libmesh/mesh.h"
65 #include "libmesh/mesh_generation.h"
66 #include "libmesh/mesh_modification.h"
67 #include "libmesh/mesh_refinement.h"
68 #include "libmesh/newton_solver.h"
69 #include "libmesh/numeric_vector.h"
70 #include "libmesh/steady_solver.h"
71 #include "libmesh/system_norm.h"
72 #include "libmesh/petsc_vector.h"
73 
74 // Adjoint Related includes
75 #include "libmesh/qoi_set.h"
76 #include "libmesh/adjoint_refinement_estimator.h"
77 
78 // libMesh I/O includes
79 #include "libmesh/getpot.h"
80 #include "libmesh/gmv_io.h"
81 
82 // Local includes
83 #include "femparameters.h"
84 #include "poisson.h"
85 
86 // Bring in everything from the libMesh namespace
87 using namespace libMesh;
88 
89 // Local function declarations
90 
91 // Number output files, the files are give a prefix of primal or adjoint_i depending on
92 // whether the output is the primal solution or the dual solution for the ith QoI
93 
94 // Write gmv output
96  unsigned int a_step, // The adaptive step count
97  std::string solution_type) // primal or adjoint solve
98 {
99 #ifdef LIBMESH_HAVE_GMV
100  MeshBase & mesh = es.get_mesh();
101 
102  std::ostringstream file_name_gmv;
103  file_name_gmv << solution_type
104  << ".out.gmv."
105  << std::setw(2)
106  << std::setfill('0')
107  << std::right
108  << a_step;
109 
111  (file_name_gmv.str(), es);
112 #endif
113 }
114 
115 // Set the parameters for the nonlinear and linear solvers to be used during the simulation
117 {
118  // Use analytical jacobians?
119  system.analytic_jacobians() = param.analytic_jacobians;
120 
121  // Verify analytic jacobians against numerical ones?
123 
124  // Use the prescribed FE type
125  system.fe_family() = param.fe_family[0];
126  system.fe_order() = param.fe_order[0];
127 
128  // More desperate debugging options
130  system.print_solutions = param.print_solutions;
132  system.print_residuals = param.print_residuals;
134  system.print_jacobians = param.print_jacobians;
135 
136  // No transient time solver
137  system.time_solver =
138  UniquePtr<TimeSolver>(new SteadySolver(system));
139 
140  // Nonlinear solver options
141  {
142  NewtonSolver * solver = new NewtonSolver(system);
143  system.time_solver->diff_solver() = UniquePtr<DiffSolver>(solver);
144 
145  solver->quiet = param.solver_quiet;
147  solver->minsteplength = param.min_step_length;
152  if (system.time_solver->reduce_deltat_on_diffsolver_failure)
153  {
154  solver->continue_after_max_iterations = true;
155  solver->continue_after_backtrack_failure = true;
156  }
157 
158  // And the linear solver options
162  }
163 }
164 
165 
166 #ifdef LIBMESH_ENABLE_AMR
167 
168 // Build the mesh refinement object and set parameters for refining/coarsening etc
170  FEMParameters & param)
171 {
172  MeshRefinement * mesh_refinement = new MeshRefinement(mesh);
173  mesh_refinement->coarsen_by_parents() = true;
174  mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
175  mesh_refinement->nelem_target() = param.nelem_target;
176  mesh_refinement->refine_fraction() = param.refine_fraction;
177  mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
178  mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
179 
180  return UniquePtr<MeshRefinement>(mesh_refinement);
181 }
182 
183 
184 // This is where declare the adjoint refined error estimator. This
185 // estimator builds an error bound for Q(u) - Q(u_h), by solving the
186 // adjoint problem on a finer Finite Element space. For more details
187 // see the description of the Adjoint Refinement Error Estimator in
188 // adjoint_refinement_error_estimator.C
190 {
191  libMesh::out << "Computing the error estimate using the Adjoint Refinement Error Estimator\n" << std::endl;
192 
193  AdjointRefinementEstimator * adjoint_refinement_estimator = new AdjointRefinementEstimator;
194 
195  adjoint_refinement_estimator->qoi_set() = qois;
196 
197  // We enrich the FE space for the dual problem by doing 2 uniform h refinements
198  adjoint_refinement_estimator->number_h_refinements = 2;
199 
200  return UniquePtr<AdjointRefinementEstimator>(adjoint_refinement_estimator);
201 }
202 
203 #endif // LIBMESH_ENABLE_AMR
204 
205 
206 // The main program.
207 int main (int argc, char ** argv)
208 {
209  // Initialize libMesh.
210  LibMeshInit init (argc, argv);
211 
212  // Skip adaptive examples on a non-adaptive libMesh build
213 #ifndef LIBMESH_ENABLE_AMR
214  libmesh_example_requires(false, "--enable-amr");
215 #else
216 
217  // This doesn't converge with Eigen BICGSTAB for some reason...
218  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc");
219 
220  libMesh::out << "Started " << argv[0] << std::endl;
221 
222  // Make sure the general input file exists, and parse it
223  {
224  std::ifstream i("general.in");
225  if (!i)
226  libmesh_error_msg('[' << init.comm().rank() << "] Can't find general.in; exiting early.");
227  }
228  GetPot infile("general.in");
229 
230  // Read in parameters from the input file
231  FEMParameters param(init.comm());
232  param.read(infile);
233 
234  // Skip this default-2D example if libMesh was compiled as 1D-only.
235  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
236 
237  // Create a mesh, with dimension to be overridden later, distributed
238  // across the default MPI communicator.
239  Mesh mesh(init.comm());
240 
241  // And an object to refine it
242  UniquePtr<MeshRefinement> mesh_refinement =
243  build_mesh_refinement(mesh, param);
244 
245  // And an EquationSystems to run on it
246  EquationSystems equation_systems (mesh);
247 
248  libMesh::out << "Reading in and building the mesh" << std::endl;
249 
250  // Read in the mesh
252  (mesh, 2, 2,
253  0.0, 1.0,
254  0.0, 1.0,
255  QUAD4);
256 
257  // Make all the elements of the mesh second order so we can compute
258  // with a higher order basis
260 
261  // Create a mesh refinement object to do the initial uniform refinements
262  // on the coarse grid read in from lshaped.xda
263  MeshRefinement initial_uniform_refinements(mesh);
264  initial_uniform_refinements.uniformly_refine(param.coarserefinements);
265 
266  libMesh::out << "Building system" << std::endl;
267 
268  // Build the FEMSystem
269  PoissonSystem & system = equation_systems.add_system<PoissonSystem> ("PoissonSystem");
270 
271  // Set its parameters
272  set_system_parameters(system, param);
273 
274  libMesh::out << "Initializing systems" << std::endl;
275 
276  equation_systems.init ();
277 
278  // Add an adjoint_solution0 vector to the system
279  system.add_vector("adjoint_solution0", false, GHOSTED);
280 
281  // Print information about the mesh and system to the screen.
282  mesh.print_info();
283  equation_systems.print_info();
284 
285  // Get a pointer to the linear solver object to be able to reuse preconditioner
286  LinearSolver<Number> * linear_solver = system.get_linear_solver();
287 
288  {
289  // Adaptively solve the timestep
290  unsigned int a_step = 0;
291  for (; a_step != param.max_adaptivesteps; ++a_step)
292  {
293  // We can't adapt to both a tolerance and a
294  // target mesh size
295  if (param.global_tolerance != 0.)
296  libmesh_assert_equal_to (param.nelem_target, 0);
297  // If we aren't adapting to a tolerance we need a
298  // target mesh size
299  else
300  libmesh_assert_greater (param.nelem_target, 0);
301 
302  // Dont reuse preconditioners before the primal solve
303  linear_solver->reuse_preconditioner(false);
304 
305  // Solve the forward problem
306  system.solve();
307 
308  // Write out the computed primal solution
309  write_output(equation_systems, a_step, "primal");
310 
311  // Declare a QoISet object, we need this object to set weights for our QoI error contributions
312  QoISet qois;
313 
314  // Declare a qoi_indices vector, each index will correspond to a QoI
315  std::vector<unsigned int> qoi_indices;
316  qoi_indices.push_back(0);
317  qois.add_indices(qoi_indices);
318 
319  // Set weights for each index, these will weight the contribution of each QoI in the final error
320  // estimate to be used for flagging elements for refinement
321  qois.set_weight(0, 1.0);
322 
323  // Make sure we get the contributions to the adjoint RHS from the sides
324  system.assemble_qoi_sides = true;
325 
326  // We are about to solve the adjoint system, but before we do this we see the same preconditioner
327  // flag to reuse the preconditioner from the forward solver
328  linear_solver->reuse_preconditioner(param.reuse_preconditioner);
329 
330  // Solve the adjoint system. This takes the transpose of the stiffness matrix and then
331  // solves the resulting system
332  system.adjoint_solve();
333 
334  //Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
335  system.set_adjoint_already_solved(true);
336 
337  // Get a pointer to the primal solution vector
338  NumericVector<Number> & primal_solution = *system.solution;
339 
340  //Get a pointer to the solution vector of the adjoint problem for QoI 0
341  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
342 
343  //Swap the primal and dual solutions so we can write out the adjoint solution
344  primal_solution.swap(dual_solution_0);
345  write_output(equation_systems, a_step, "adjoint_0");
346 
347  //Swap back
348  primal_solution.swap(dual_solution_0);
349 
350  libMesh::out << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
351  << " active elements and "
352  << equation_systems.n_active_dofs()
353  << " active dofs." << std::endl ;
354 
355  // Postprocess, compute the approximate QoIs and write them out to the console
356  libMesh::out << "Postprocessing: " << std::endl;
357  system.postprocess_sides = true;
358  system.postprocess();
359 
360  Number QoI_0_computed = system.get_QoI_value("computed", 0);
361  Number QoI_0_exact = system.get_QoI_value("exact", 0);
362 
363  libMesh::out << "The computed QoI 0 is " << std::setprecision(17)
364  << QoI_0_computed << std::endl;
365  libMesh::out << "The relative error in QoI 0 is " << std::setprecision(17)
366  << std::abs(QoI_0_computed - QoI_0_exact) << std::endl; // / std::abs(QoI_0_exact)
367 
368  // We will declare an error vector for passing to the adjoint refinement error estimator
369  ErrorVector QoI_elementwise_error;
370 
371  // Build an adjoint refinement error estimator object
372  UniquePtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
374 
375  // Estimate the error in each element using the Adjoint Refinement estimator
376  adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
377 
378  // Print out the computed error estimate, note that we access the global error estimates
379  // using an accessor function, right now sum(QoI_elementwise_error) != global_QoI_error_estimate
380  libMesh::out << "The computed relative error in QoI 0 is " << std::setprecision(17)
381  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) << std::endl; // / std::abs(QoI_0_exact)
382 
383  // Also print out effectivity indices (estimated error/true error)
384  libMesh::out << "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
385  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
386  std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
387 
388  // For refinement purposes we need to sort by error
389  // *magnitudes*, but AdjointRefinement gives us signed errors.
390  if (!param.refine_uniformly)
391  for (std::size_t i=0; i<QoI_elementwise_error.size(); i++)
392  if (QoI_elementwise_error[i] != 0.)
393  QoI_elementwise_error[i] = std::abs(QoI_elementwise_error[i]);
394 
395  // We have to refine either based on reaching an error tolerance or
396  // a number of elements target, which should be verified above
397  // Otherwise we flag elements by error tolerance or nelem target
398 
399  // Uniform refinement
400  if (param.refine_uniformly)
401  {
402  mesh_refinement->uniformly_refine(1);
403  }
404  // Adaptively refine based on reaching an error tolerance
405  else if (param.global_tolerance >= 0. && param.nelem_target == 0.)
406  {
407  mesh_refinement->flag_elements_by_error_tolerance (QoI_elementwise_error);
408 
409  mesh_refinement->refine_and_coarsen_elements();
410  }
411  // Adaptively refine based on reaching a target number of elements
412  else
413  {
414  if (mesh.n_active_elem() >= param.nelem_target)
415  {
416  libMesh::out << "We reached the target number of elements.\n" << std::endl;
417  break;
418  }
419 
420  mesh_refinement->flag_elements_by_nelem_target (QoI_elementwise_error);
421 
422  mesh_refinement->refine_and_coarsen_elements();
423  }
424 
425  // Dont forget to reinit the system after each adaptive refinement!
426  equation_systems.reinit();
427 
428  libMesh::out << "Refined mesh to "
429  << mesh.n_active_elem()
430  << " active elements and "
431  << equation_systems.n_active_dofs()
432  << " active dofs." << std::endl;
433  }
434 
435  // On the last adaptive step, dont refine elements and check regressions via asserts
436  if (a_step == param.max_adaptivesteps)
437  {
438  linear_solver->reuse_preconditioner(false);
439  system.solve();
440 
441  write_output(equation_systems, a_step, "primal");
442 
443  NumericVector<Number> & primal_solution = *system.solution;
444 
445  QoISet qois;
446  std::vector<unsigned int> qoi_indices;
447 
448  qoi_indices.push_back(0);
449  qois.add_indices(qoi_indices);
450 
451  qois.set_weight(0, 1.0);
452 
453  system.assemble_qoi_sides = true;
454 
455  linear_solver->reuse_preconditioner(param.reuse_preconditioner);
456  system.adjoint_solve();
457  system.set_adjoint_already_solved(true);
458 
459  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
460 
461  primal_solution.swap(dual_solution_0);
462  write_output(equation_systems, a_step, "adjoint_0");
463 
464  primal_solution.swap(dual_solution_0);
465 
466  libMesh::out << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
467  << " active elements and "
468  << equation_systems.n_active_dofs()
469  << " active dofs." << std::endl ;
470 
471  libMesh::out << "Postprocessing: " << std::endl;
472  system.postprocess_sides = true;
473  system.postprocess();
474 
475  Number QoI_0_computed = system.get_QoI_value("computed", 0);
476  Number QoI_0_exact = system.get_QoI_value("exact", 0);
477 
478  libMesh::out << "The computed QoI 0 is " << std::setprecision(17)
479  << QoI_0_computed << std::endl;
480  libMesh::out << "The relative error in QoI 0 is " << std::setprecision(17)
481  << std::abs(QoI_0_computed - QoI_0_exact) << std::endl; // / std::abs(QoI_0_exact)
482 
483 
484  // We will declare an error vector for passing to the adjoint
485  // refinement error estimator Right now, only the first entry
486  // of this vector will be filled (with the global QoI error
487  // estimate) Later, each entry of the vector will contain
488  // elementwise error that the user can sum to get the total
489  // error
490  ErrorVector QoI_elementwise_error;
491 
492  // Build an adjoint refinement error estimator object
493  UniquePtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
495 
496  // Estimate the error in each element using the Adjoint Refinement estimator
497  adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
498 
499  // Print out the computed error estimate, note that we access
500  // the global error estimates using an accessor function,
501  // right now sum(QoI_elementwise_error) != global_QoI_error_estimate
502  libMesh::out << "The computed relative error in QoI 0 is " << std::setprecision(17)
503  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) << std::endl; // / std::abs(QoI_0_exact)
504 
505  // Also print out effectivity indices (estimated error/true error)
506  libMesh::out << "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
507  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
508  std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
509 
510  // Hard coded assert to ensure that the actual numbers we are getting are what they should be
511 
512  // The effectivity index isn't exactly reproducible at single precision
513  // libmesh_assert_less(std::abs(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact) - 0.84010976704434637), 1.e-5);
514  // libmesh_assert_less(std::abs(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact) - 0.48294428289950514), 1.e-5);
515 
516  // But the effectivity indices should always be sane
517  // libmesh_assert_less(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact), 2.5);
518  // libmesh_assert_greater(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact), .4);
519  // libmesh_assert_less(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact), 2.5);
520  // libmesh_assert_greater(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact), .4);
521 
522  // And the computed errors should still be low
523  // libmesh_assert_less(std::abs(QoI_0_computed - QoI_0_exact), 2e-4);
524  // libmesh_assert_less(std::abs(QoI_1_computed - QoI_1_exact), 2e-4);
525  }
526  }
527 
528  libMesh::err << '[' << mesh.processor_id()
529  << "] Completing output." << std::endl;
530 
531 #endif // #ifndef LIBMESH_ENABLE_AMR
532 
533  // All done.
534  return 0;
535 }
bool continue_after_max_iterations
Defaults to true, telling the DiffSolver to continue rather than exit when a solve has reached its ma...
Definition: diff_solver.h:174
unsigned int nelem_target
Definition: femparameters.h:56
OStreamProxy err
bool print_solution_norms
double abs(double a)
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
Real minimum_linear_tolerance
The tolerance for linear solves is kept above this minimum.
Definition: diff_solver.h:215
QoISet & qoi_set()
Access to the QoISet (default: weight all QoIs equally) to use when computing errors.
This is the EquationSystems class.
libMesh::Real initial_linear_tolerance
virtual dof_id_type n_active_elem() const =0
Real verify_analytic_jacobians
If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calc...
Definition: fem_system.h:203
void set_weight(unsigned int, Real)
Set the weight for this index.
Definition: qoi_set.h:229
void set_adjoint_already_solved(bool setting)
Setter for the adjoint_already_solved boolean.
Definition: system.h:378
bool analytic_jacobians
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:162
This class implements a ``brute force&#39;&#39; goal-oriented error estimator which computes an estimate of e...
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Definition: diff_system.h:329
Real & absolute_global_tolerance()
If absolute_global_tolerance is set to a nonzero value, methods like flag_elements_by_global_toleranc...
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) libmesh_override
This function sets the _is_adjoint boolean member of TimeSolver to true and then calls the adjoint_so...
Definition: diff_system.C:164
unsigned int max_nonlinear_iterations
The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_...
Definition: diff_solver.h:156
bool require_residual_reduction
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
virtual void postprocess(void)
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: poisson.C:186
MeshBase & mesh
bool print_residual_norms
Real & refine_fraction()
The refine_fraction sets either a desired target or a desired maximum number of elements to flag for ...
libMesh::Real relative_step_tolerance
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
bool postprocess_sides
If postprocess_sides is true (it is false by default), the postprocessing loop will loop over all sid...
Definition: diff_system.h:302
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.
Real initial_linear_tolerance
Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the ...
Definition: diff_solver.h:210
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
libMesh::Real refine_fraction
Definition: femparameters.h:58
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
Definition: diff_system.C:175
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:334
bool & coarsen_by_parents()
If coarsen_by_parents is true, complete groups of sibling elements (elements with the same parent) wi...
This class defines a solver which uses the default libMesh linear solver in a quasiNewton method to h...
Definition: newton_solver.h:46
This is the MeshBase class.
Definition: mesh_base.h:68
std::vector< unsigned int > fe_order
UniquePtr< MeshRefinement > build_mesh_refinement(MeshBase &mesh, FEMParameters &param)
Definition: adjoints_ex6.C:169
dof_id_type & nelem_target()
If nelem_target is set to a nonzero value, methods like flag_elements_by_nelem_target() will attempt ...
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscDiffSolver & solver
SolverPackage default_solver_package()
Definition: libmesh.C:995
libMesh::Real coarsen_threshold
Definition: femparameters.h:58
bool print_solution_norms
Set print_residual_norms to true to print |U| whenever it is used in an assembly() call...
Definition: diff_system.h:308
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
This is the MeshRefinement class.
unsigned int max_linear_iterations
bool print_residual_norms
Set print_residual_norms to true to print |F| whenever it is assembled.
Definition: diff_system.h:319
virtual void all_second_order(const bool full_ordered=true)=0
Converts a (conforming, non-refined) mesh with linear elements into a mesh with second-order elements...
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:148
libMesh::Real minimum_linear_tolerance
virtual void reinit()
Reinitialize all the systems.
std::vector< std::string > fe_family
void write_output(EquationSystems &es, unsigned int a_step, std::string solution_type)
Definition: adjoints_ex6.C:95
libMesh::Real verify_analytic_jacobians
Real & coarsen_threshold()
The coarsen_threshold provides hysteresis in AMR/C strategies.
bool print_residuals
Set print_residuals to true to print F whenever it is assembled.
Definition: diff_system.h:324
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
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.
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
libMesh::Real global_tolerance
Definition: femparameters.h:57
Real linear_tolerance_multiplier
The tolerance for linear solves is kept below this multiplier (which defaults to 1e-3) times the norm...
std::string & fe_family()
Definition: poisson.h:21
bool print_solutions
Set print_solutions to true to print U whenever it is used in an assembly() call. ...
Definition: diff_system.h:314
libMesh::Real relative_residual_tolerance
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
unsigned int & fe_order()
Definition: poisson.h:22
bool continue_after_backtrack_failure
Defaults to false, telling the DiffSolver to throw an error when the backtracking scheme fails to fin...
Definition: diff_solver.h:180
virtual void reuse_preconditioner(bool)
Set the same_preconditioner flag, which indicates if we reuse the same preconditioner for subsequent ...
bool print_jacobian_norms
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
Real minsteplength
If the quasi-Newton step length must be reduced to below this factor to give a residual reduction...
bool assemble_qoi_sides
If assemble_qoi_sides is true (it is false by default), the assembly loop for a quantity of interest ...
Definition: diff_qoi.h:82
OStreamProxy out
void add_indices(const std::vector< unsigned int > &indices)
Add this indices to the set to be calculated.
Definition: qoi_set.C:46
libMesh::Real min_step_length
void read(GetPot &input, const std::vector< std::string > *other_variable_names=libmesh_nullptr)
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:989
const MeshBase & get_mesh() const
unsigned int max_nonlinear_iterations
UniquePtr< AdjointRefinementEstimator > build_adjoint_refinement_error_estimator(QoISet &qois)
Definition: adjoints_ex6.C:189
unsigned int rank() const
Definition: parallel.h:724
bool require_residual_reduction
If this is set to true, the solver is forced to test the residual after each Newton step...
Definition: newton_solver.h:98
bool & analytic_jacobians()
Definition: poisson.h:23
virtual void init()
Initialize all the systems.
void set_system_parameters(PoissonSystem &system, FEMParameters &param)
Definition: adjoints_ex6.C:116
std::size_t n_active_dofs() const
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
int main(int argc, char **argv)
Definition: adjoints_ex6.C:207
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
libMesh::Real linear_tolerance_multiplier
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
Number & get_QoI_value(std::string type, unsigned int QoI_index)
Definition: poisson.h:29
Real relative_residual_tolerance
Definition: diff_solver.h:192
libMesh::Real coarsen_fraction
Definition: femparameters.h:58
processor_id_type processor_id() const
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.