libMesh
Functions
adjoints_ex4.C File Reference

Go to the source code of this file.

Functions

void write_output (EquationSystems &es, unsigned int a_step, std::string solution_type, FEMParameters &param)
 
void set_system_parameters (LaplaceSystem &system, FEMParameters &param)
 
UniquePtr< MeshRefinementbuild_mesh_refinement (MeshBase &mesh, FEMParameters &param)
 
UniquePtr< AdjointRefinementEstimatorbuild_adjoint_refinement_error_estimator (QoISet &qois)
 
int main (int argc, char **argv)
 

Function Documentation

UniquePtr<AdjointRefinementEstimator> build_adjoint_refinement_error_estimator ( QoISet qois)

Definition at line 227 of file adjoints_ex4.C.

References libMesh::AdjointRefinementEstimator::number_h_refinements, libMesh::out, and libMesh::AdjointRefinementEstimator::qoi_set().

Referenced by main().

228 {
229  libMesh::out << "Computing the error estimate using the Adjoint Refinement Error Estimator" << std::endl << std::endl;
230 
231  AdjointRefinementEstimator *adjoint_refinement_estimator = new AdjointRefinementEstimator;
232 
233  adjoint_refinement_estimator->qoi_set() = qois;
234 
235  // We enrich the FE space for the dual problem by doing 2 uniform h refinements
236  adjoint_refinement_estimator->number_h_refinements = 2;
237 
238  return UniquePtr<AdjointRefinementEstimator>(adjoint_refinement_estimator);
239 }
QoISet & qoi_set()
Access to the QoISet (default: weight all QoIs equally) to use when computing errors.
This class implements a ``brute force&#39;&#39; goal-oriented error estimator which computes an estimate of e...
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
OStreamProxy out
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
UniquePtr<MeshRefinement> build_mesh_refinement ( MeshBase mesh,
FEMParameters param 
)

Definition at line 210 of file adjoints_ex4.C.

References libMesh::MeshRefinement::absolute_global_tolerance(), libMesh::MeshRefinement::coarsen_by_parents(), FEMParameters::coarsen_fraction, libMesh::MeshRefinement::coarsen_fraction(), FEMParameters::coarsen_threshold, libMesh::MeshRefinement::coarsen_threshold(), FEMParameters::global_tolerance, FEMParameters::nelem_target, libMesh::MeshRefinement::nelem_target(), FEMParameters::refine_fraction, and libMesh::MeshRefinement::refine_fraction().

Referenced by main().

212 {
213  MeshRefinement * mesh_refinement = new MeshRefinement(mesh);
214  mesh_refinement->coarsen_by_parents() = true;
215  mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
216  mesh_refinement->nelem_target() = param.nelem_target;
217  mesh_refinement->refine_fraction() = param.refine_fraction;
218  mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
219  mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
220 
221  return UniquePtr<MeshRefinement>(mesh_refinement);
222 }
unsigned int nelem_target
Definition: femparameters.h:56
Real & absolute_global_tolerance()
If absolute_global_tolerance is set to a nonzero value, methods like flag_elements_by_global_toleranc...
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
Real & refine_fraction()
The refine_fraction sets either a desired target or a desired maximum number of elements to flag for ...
libMesh::Real refine_fraction
Definition: femparameters.h:58
bool & coarsen_by_parents()
If coarsen_by_parents is true, complete groups of sibling elements (elements with the same parent) wi...
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
libMesh::Real coarsen_threshold
Definition: femparameters.h:58
This is the MeshRefinement class.
Real & coarsen_threshold()
The coarsen_threshold provides hysteresis in AMR/C strategies.
libMesh::Real global_tolerance
Definition: femparameters.h:57
libMesh::Real coarsen_fraction
Definition: femparameters.h:58
int main ( int  argc,
char **  argv 
)

Definition at line 245 of file adjoints_ex4.C.

References std::abs(), libMesh::QoISet::add_indices(), libMesh::EquationSystems::add_system(), libMesh::DifferentiableSystem::adjoint_solve(), libMesh::MeshBase::all_second_order(), libMesh::DifferentiableQoI::assemble_qoi_sides, build_adjoint_refinement_error_estimator(), build_mesh_refinement(), libMesh::LibMeshInit::comm(), libMesh::default_solver_package(), libMesh::EIGEN_SOLVERS, libMesh::err, libMesh::System::get_adjoint_solution(), libMesh::DifferentiableSystem::get_linear_solver(), LaplaceSystem::get_QoI_value(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::EquationSystems::n_active_dofs(), libMesh::MeshBase::n_active_elem(), libMesh::out, LaplaceSystem::postprocess(), libMesh::DifferentiableSystem::postprocess_sides, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::ParallelObject::processor_id(), libMesh::Parallel::Communicator::rank(), FEMParameters::read(), libMesh::MeshBase::read(), libMesh::EquationSystems::reinit(), libMesh::LinearSolver< T >::reuse_preconditioner(), libMesh::System::set_adjoint_already_solved(), set_system_parameters(), libMesh::QoISet::set_weight(), libMesh::System::solution, libMesh::FEMSystem::solve(), libMesh::NumericVector< T >::swap(), libMesh::MeshRefinement::uniformly_refine(), and write_output().

246 {
247  // Initialize libMesh.
248  LibMeshInit init (argc, argv);
249 
250  // This example requires a linear solver package.
251  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
252  "--enable-petsc, --enable-trilinos, or --enable-eigen");
253 
254  // Skip adaptive examples on a non-adaptive libMesh build
255 #ifndef LIBMESH_ENABLE_AMR
256  libmesh_example_requires(false, "--enable-amr");
257 #else
258 
259  // This doesn't converge with Eigen BICGSTAB for some reason...
260  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc");
261 
262  libMesh::out << "Started " << argv[0] << std::endl;
263 
264  // Make sure the general input file exists, and parse it
265  {
266  std::ifstream i("general.in");
267  if (!i)
268  libmesh_error_msg('[' << init.comm().rank() << "] Can't find general.in; exiting early.");
269  }
270  GetPot infile("general.in");
271 
272  // Read in parameters from the input file
273  FEMParameters param(init.comm());
274  param.read(infile);
275 
276  // Skip this default-2D example if libMesh was compiled as 1D-only.
277  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
278 
279  // Create a mesh, with dimension to be overridden later, distributed
280  // across the default MPI communicator.
281  Mesh mesh(init.comm());
282 
283  // And an object to refine it
284  UniquePtr<MeshRefinement> mesh_refinement =
285  build_mesh_refinement(mesh, param);
286 
287  // And an EquationSystems to run on it
288  EquationSystems equation_systems (mesh);
289 
290  libMesh::out << "Reading in and building the mesh" << std::endl;
291 
292  // Read in the mesh
293  mesh.read(param.domainfile.c_str());
294  // Make all the elements of the mesh second order so we can compute
295  // with a higher order basis
297 
298  // Create a mesh refinement object to do the initial uniform refinements
299  // on the coarse grid read in from lshaped.xda
300  MeshRefinement initial_uniform_refinements(mesh);
301  initial_uniform_refinements.uniformly_refine(param.coarserefinements);
302 
303  libMesh::out << "Building system" << std::endl;
304 
305  // Build the FEMSystem
306  LaplaceSystem & system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
307 
308  // Set its parameters
309  set_system_parameters(system, param);
310 
311  libMesh::out << "Initializing systems" << std::endl;
312 
313  equation_systems.init ();
314 
315  // Print information about the mesh and system to the screen.
316  mesh.print_info();
317  equation_systems.print_info();
318  LinearSolver<Number> *linear_solver = system.get_linear_solver();
319 
320  {
321  // Adaptively solve the timestep
322  unsigned int a_step = 0;
323  for (; a_step != param.max_adaptivesteps; ++a_step)
324  {
325  // We can't adapt to both a tolerance and a
326  // target mesh size
327  if (param.global_tolerance != 0.)
328  libmesh_assert_equal_to (param.nelem_target, 0);
329  // If we aren't adapting to a tolerance we need a
330  // target mesh size
331  else
332  libmesh_assert_greater (param.nelem_target, 0);
333 
334  linear_solver->reuse_preconditioner(false);
335 
336  // Solve the forward problem
337  system.solve();
338 
339  // Write out the computed primal solution
340  write_output(equation_systems, a_step, "primal", param);
341 
342  // Get a pointer to the primal solution vector
343  NumericVector<Number> & primal_solution = *system.solution;
344 
345  // Declare a QoISet object, we need this object to set weights for our QoI error contributions
346  QoISet qois;
347 
348  // Declare a qoi_indices vector, each index will correspond to a QoI
349  std::vector<unsigned int> qoi_indices;
350  qoi_indices.push_back(0);
351  qoi_indices.push_back(1);
352  qois.add_indices(qoi_indices);
353 
354  // Set weights for each index, these will weight the contribution of each QoI in the final error
355  // estimate to be used for flagging elements for refinement
356  qois.set_weight(0, 0.5);
357  qois.set_weight(1, 0.5);
358 
359  // Make sure we get the contributions to the adjoint RHS from the sides
360  system.assemble_qoi_sides = true;
361 
362  // We are about to solve the adjoint system, but before we do this we see the same preconditioner
363  // flag to reuse the preconditioner from the forward solver
364  linear_solver->reuse_preconditioner(param.reuse_preconditioner);
365 
366  // Solve the adjoint system. This takes the transpose of the stiffness matrix and then
367  // solves the resulting system
368  system.adjoint_solve();
369 
370  // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
371  system.set_adjoint_already_solved(true);
372 
373  // Get a pointer to the solution vector of the adjoint problem for QoI 0
374  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
375 
376  // Swap the primal and dual solutions so we can write out the adjoint solution
377  primal_solution.swap(dual_solution_0);
378  write_output(equation_systems, a_step, "adjoint_0", param);
379 
380  // Swap back
381  primal_solution.swap(dual_solution_0);
382 
383  // Get a pointer to the solution vector of the adjoint problem for QoI 0
384  NumericVector<Number> & dual_solution_1 = system.get_adjoint_solution(1);
385 
386  // Swap again
387  primal_solution.swap(dual_solution_1);
388  write_output(equation_systems, a_step, "adjoint_1", param);
389 
390  // Swap back again
391  primal_solution.swap(dual_solution_1);
392 
393  libMesh::out << "Adaptive step "
394  << a_step
395  << ", we have "
396  << mesh.n_active_elem()
397  << " active elements and "
398  << equation_systems.n_active_dofs()
399  << " active dofs."
400  << std::endl;
401 
402  // Postprocess, compute the approximate QoIs and write them out to the console
403  libMesh::out << "Postprocessing: " << std::endl;
404  system.postprocess_sides = true;
405  system.postprocess();
406  Number QoI_0_computed = system.get_QoI_value("computed", 0);
407  Number QoI_0_exact = system.get_QoI_value("exact", 0);
408  Number QoI_1_computed = system.get_QoI_value("computed", 1);
409  Number QoI_1_exact = system.get_QoI_value("exact", 1);
410 
411  libMesh::out << "The relative error in QoI 0 is "
412  << std::setprecision(17)
413  << std::abs(QoI_0_computed - QoI_0_exact) / std::abs(QoI_0_exact)
414  << std::endl;
415 
416  libMesh::out << "The relative error in QoI 1 is "
417  << std::setprecision(17)
418  << std::abs(QoI_1_computed - QoI_1_exact) / std::abs(QoI_1_exact)
419  << std::endl
420  << std::endl;
421 
422  // We will declare an error vector for passing to the adjoint refinement error estimator
423  ErrorVector QoI_elementwise_error;
424 
425  // Build an adjoint refinement error estimator object
426  UniquePtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
428 
429  // Estimate the error in each element using the Adjoint Refinement estimator
430  adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
431 
432  // Print out the computed error estimate, note that we access the global error estimates
433  // using an accessor function, right now sum(QoI_elementwise_error) != global_QoI_error_estimate
434  libMesh::out << "The computed relative error in QoI 0 is "
435  << std::setprecision(17)
436  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_exact)
437  << std::endl;
438 
439  libMesh::out << "The computed relative error in QoI 1 is "
440  << std::setprecision(17)
441  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_exact)
442  << std::endl
443  << std::endl;
444 
445  // Also print out effectivity indices (estimated error/true error)
446  libMesh::out << "The effectivity index for the computed error in QoI 0 is "
447  << std::setprecision(17)
448  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact)
449  << std::endl;
450 
451  libMesh::out << "The effectivity index for the computed error in QoI 1 is "
452  << std::setprecision(17)
453  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact)
454  << std::endl
455  << std::endl;
456 
457  // For refinement purposes we need to sort by error
458  // *magnitudes*, but AdjointRefinement gives us signed errors.
459  if (!param.refine_uniformly)
460  for (std::size_t i=0; i<QoI_elementwise_error.size(); i++)
461  if (QoI_elementwise_error[i] != 0.)
462  QoI_elementwise_error[i] = std::abs(QoI_elementwise_error[i]);
463 
464  // We have to refine either based on reaching an error tolerance or
465  // a number of elements target, which should be verified above
466  // Otherwise we flag elements by error tolerance or nelem target
467 
468  // Uniform refinement
469  if (param.refine_uniformly)
470  {
471  mesh_refinement->uniformly_refine(1);
472  }
473  // Adaptively refine based on reaching an error tolerance
474  else if (param.global_tolerance >= 0. && param.nelem_target == 0.)
475  {
476  mesh_refinement->flag_elements_by_error_tolerance (QoI_elementwise_error);
477 
478  mesh_refinement->refine_and_coarsen_elements();
479  }
480  // Adaptively refine based on reaching a target number of elements
481  else
482  {
483  if (mesh.n_active_elem() >= param.nelem_target)
484  {
485  libMesh::out << "We reached the target number of elements." << std::endl << std::endl;
486  break;
487  }
488 
489  mesh_refinement->flag_elements_by_nelem_target (QoI_elementwise_error);
490 
491  mesh_refinement->refine_and_coarsen_elements();
492  }
493 
494  // Dont forget to reinit the system after each adaptive refinement !
495  equation_systems.reinit();
496 
497  libMesh::out << "Refined mesh to "
498  << mesh.n_active_elem()
499  << " active elements and "
500  << equation_systems.n_active_dofs()
501  << " active dofs."
502  << std::endl;
503  }
504 
505  // Do one last solve if necessary
506  if (a_step == param.max_adaptivesteps)
507  {
508  linear_solver->reuse_preconditioner(false);
509  system.solve();
510 
511  write_output(equation_systems, a_step, "primal", param);
512 
513  NumericVector<Number> & primal_solution = *system.solution;
514 
515  QoISet qois;
516  std::vector<unsigned int> qoi_indices;
517 
518  qoi_indices.push_back(0);
519  qoi_indices.push_back(1);
520  qois.add_indices(qoi_indices);
521 
522  qois.set_weight(0, 0.5);
523  qois.set_weight(1, 0.5);
524 
525  system.assemble_qoi_sides = true;
526  linear_solver->reuse_preconditioner(param.reuse_preconditioner);
527  system.adjoint_solve();
528 
529  // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
530  system.set_adjoint_already_solved(true);
531 
532  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
533 
534  primal_solution.swap(dual_solution_0);
535  write_output(equation_systems, a_step, "adjoint_0", param);
536 
537  primal_solution.swap(dual_solution_0);
538 
539  NumericVector<Number> & dual_solution_1 = system.get_adjoint_solution(1);
540 
541  primal_solution.swap(dual_solution_1);
542  write_output(equation_systems, a_step, "adjoint_1", param);
543 
544  primal_solution.swap(dual_solution_1);
545 
546  libMesh::out << "Adaptive step "
547  << a_step
548  << ", we have "
549  << mesh.n_active_elem()
550  << " active elements and "
551  << equation_systems.n_active_dofs()
552  << " active dofs."
553  << std::endl;
554 
555  libMesh::out << "Postprocessing: " << std::endl;
556  system.postprocess_sides = true;
557  system.postprocess();
558 
559  Number QoI_0_computed = system.get_QoI_value("computed", 0);
560  Number QoI_0_exact = system.get_QoI_value("exact", 0);
561  Number QoI_1_computed = system.get_QoI_value("computed", 1);
562  Number QoI_1_exact = system.get_QoI_value("exact", 1);
563 
564  libMesh::out << "The relative error in QoI 0 is "
565  << std::setprecision(17)
566  << std::abs(QoI_0_computed - QoI_0_exact) / std::abs(QoI_0_exact)
567  << std::endl;
568 
569  libMesh::out << "The relative error in QoI 1 is "
570  << std::setprecision(17)
571  << std::abs(QoI_1_computed - QoI_1_exact) / std::abs(QoI_1_exact)
572  << std::endl
573  << std::endl;
574 
575  // We will declare an error vector for passing to the adjoint refinement error estimator
576  // Right now, only the first entry of this vector will be filled (with the global QoI error estimate)
577  // Later, each entry of the vector will contain elementwise error that the user can sum to get the total error
578  ErrorVector QoI_elementwise_error;
579 
580  // Build an adjoint refinement error estimator object
581  UniquePtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
583 
584  // Estimate the error in each element using the Adjoint Refinement estimator
585  adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
586 
587  // Print out the computed error estimate, note that we access the global error estimates
588  // using an accessor function, right now sum(QoI_elementwise_error) != global_QoI_error_estimate
589  libMesh::out << "The computed relative error in QoI 0 is "
590  << std::setprecision(17)
591  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_exact)
592  << std::endl;
593 
594  libMesh::out << "The computed relative error in QoI 1 is "
595  << std::setprecision(17)
596  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_exact)
597  << std::endl
598  << std::endl;
599 
600  // Also print out effectivity indices (estimated error/true error)
601  libMesh::out << "The effectivity index for the computed error in QoI 0 is "
602  << std::setprecision(17)
603  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact)
604  << std::endl;
605 
606  libMesh::out << "The effectivity index for the computed error in QoI 1 is "
607  << std::setprecision(17)
608  << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact)
609  << std::endl
610  << std::endl;
611 
612  // Hard coded assert to ensure that the actual numbers we are getting are what they should be
613 
614  // The effectivity index isn't exactly reproducible at single precision
615  // 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);
616  // 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);
617 
618  // But the effectivity indices should always be sane
619  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);
620  libmesh_assert_greater(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact), .4);
621  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);
622  libmesh_assert_greater(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact), .4);
623 
624  // And the computed errors should still be low
625  libmesh_assert_less(std::abs(QoI_0_computed - QoI_0_exact), 2e-4);
626  libmesh_assert_less(std::abs(QoI_1_computed - QoI_1_exact), 2e-4);
627  }
628  }
629 
630  libMesh::err << '[' << mesh.processor_id()
631  << "] Completing output."
632  << std::endl;
633 
634 #endif // #ifndef LIBMESH_ENABLE_AMR
635 
636  // All done.
637  return 0;
638 }
OStreamProxy err
double abs(double a)
This is the EquationSystems class.
virtual dof_id_type n_active_elem() const =0
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
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
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
Number & get_QoI_value(std::string type, unsigned int QoI_index)
Definition: L-shaped.h:32
void write_output(EquationSystems &es, unsigned int a_step, std::string solution_type, FEMParameters &param)
Definition: adjoints_ex4.C:92
MeshBase & mesh
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
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
Definition: diff_system.C:175
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
This is the MeshRefinement class.
void set_system_parameters(LaplaceSystem &system, FEMParameters &param)
Definition: adjoints_ex4.C:156
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...
UniquePtr< AdjointRefinementEstimator > build_adjoint_refinement_error_estimator(QoISet &qois)
Definition: adjoints_ex4.C:227
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void reuse_preconditioner(bool)
Set the same_preconditioner flag, which indicates if we reuse the same preconditioner for subsequent ...
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
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
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
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
UniquePtr< MeshRefinement > build_mesh_refinement(MeshBase &mesh, FEMParameters &param)
Definition: adjoints_ex4.C:210
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.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
virtual void postprocess()
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: L-shaped.C:168
processor_id_type processor_id() const
void set_system_parameters ( LaplaceSystem system,
FEMParameters param 
)

Definition at line 156 of file adjoints_ex4.C.

References LaplaceSystem::analytic_jacobians(), FEMParameters::analytic_jacobians, libMesh::DiffSolver::continue_after_backtrack_failure, libMesh::DiffSolver::continue_after_max_iterations, LaplaceSystem::fe_family(), FEMParameters::fe_family, LaplaceSystem::fe_order(), FEMParameters::fe_order, FEMParameters::initial_linear_tolerance, libMesh::DiffSolver::initial_linear_tolerance, FEMParameters::linear_tolerance_multiplier, libMesh::NewtonSolver::linear_tolerance_multiplier, FEMParameters::max_linear_iterations, libMesh::DiffSolver::max_linear_iterations, FEMParameters::max_nonlinear_iterations, libMesh::DiffSolver::max_nonlinear_iterations, FEMParameters::min_step_length, FEMParameters::minimum_linear_tolerance, libMesh::DiffSolver::minimum_linear_tolerance, libMesh::NewtonSolver::minsteplength, FEMParameters::print_jacobian_norms, libMesh::DifferentiableSystem::print_jacobian_norms, FEMParameters::print_jacobians, libMesh::DifferentiableSystem::print_jacobians, FEMParameters::print_residual_norms, libMesh::DifferentiableSystem::print_residual_norms, FEMParameters::print_residuals, libMesh::DifferentiableSystem::print_residuals, FEMParameters::print_solution_norms, libMesh::DifferentiableSystem::print_solution_norms, FEMParameters::print_solutions, libMesh::DifferentiableSystem::print_solutions, libMesh::DiffSolver::quiet, FEMParameters::relative_residual_tolerance, libMesh::DiffSolver::relative_residual_tolerance, FEMParameters::relative_step_tolerance, libMesh::DiffSolver::relative_step_tolerance, libMesh::NewtonSolver::require_residual_reduction, FEMParameters::require_residual_reduction, libMesh::solver, FEMParameters::solver_quiet, libMesh::DifferentiableSystem::time_solver, FEMParameters::verify_analytic_jacobians, and libMesh::FEMSystem::verify_analytic_jacobians.

Referenced by main().

158 {
159  // Use analytical jacobians?
160  system.analytic_jacobians() = param.analytic_jacobians;
161 
162  // Verify analytic jacobians against numerical ones?
164 
165  // Use the prescribed FE type
166  system.fe_family() = param.fe_family[0];
167  system.fe_order() = param.fe_order[0];
168 
169  // More desperate debugging options
171  system.print_solutions = param.print_solutions;
173  system.print_residuals = param.print_residuals;
175  system.print_jacobians = param.print_jacobians;
176 
177  // No transient time solver
178  system.time_solver =
179  UniquePtr<TimeSolver>(new SteadySolver(system));
180 
181  // Nonlinear solver options
182  {
183  NewtonSolver * solver = new NewtonSolver(system);
184  system.time_solver->diff_solver() = UniquePtr<DiffSolver>(solver);
185 
186  solver->quiet = param.solver_quiet;
188  solver->minsteplength = param.min_step_length;
193  if (system.time_solver->reduce_deltat_on_diffsolver_failure)
194  {
195  solver->continue_after_max_iterations = true;
196  solver->continue_after_backtrack_failure = true;
197  }
198 
199  // And the linear solver options
203  }
204 }
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
bool print_solution_norms
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
libMesh::Real initial_linear_tolerance
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
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
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Definition: diff_system.h:329
std::string & fe_family()
Definition: L-shaped.h:24
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
bool print_residual_norms
libMesh::Real relative_step_tolerance
bool & analytic_jacobians()
Definition: L-shaped.h:26
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
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:334
This class defines a solver which uses the default libMesh linear solver in a quasiNewton method to h...
Definition: newton_solver.h:46
std::vector< unsigned int > fe_order
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscDiffSolver & solver
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
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
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
std::vector< std::string > fe_family
libMesh::Real verify_analytic_jacobians
bool print_residuals
Set print_residuals to true to print F whenever it is assembled.
Definition: diff_system.h:324
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
Real linear_tolerance_multiplier
The tolerance for linear solves is kept below this multiplier (which defaults to 1e-3) times the norm...
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
unsigned int & fe_order()
Definition: L-shaped.h:25
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
bool print_jacobian_norms
Real minsteplength
If the quasi-Newton step length must be reduced to below this factor to give a residual reduction...
libMesh::Real min_step_length
unsigned int max_nonlinear_iterations
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
libMesh::Real linear_tolerance_multiplier
Real relative_residual_tolerance
Definition: diff_solver.h:192
void write_output ( EquationSystems es,
unsigned int  a_step,
std::string  solution_type,
FEMParameters param 
)

Definition at line 92 of file adjoints_ex4.C.

References libMesh::EquationSystems::get_mesh(), libMesh::libmesh_ignore(), mesh, FEMParameters::output_exodus, FEMParameters::output_gmv, libMesh::MeshOutput< MT >::write_equation_systems(), and libMesh::ExodusII_IO::write_timestep().

Referenced by main().

96 {
97  // Ignore parameters when there are no output formats available.
98  libmesh_ignore(es);
99  libmesh_ignore(a_step);
100  libmesh_ignore(solution_type);
101  libmesh_ignore(param);
102 
103 #ifdef LIBMESH_HAVE_GMV
104  if (param.output_gmv)
105  {
106  MeshBase & mesh = es.get_mesh();
107 
108  std::ostringstream file_name_gmv;
109  file_name_gmv << solution_type
110  << ".out.gmv."
111  << std::setw(2)
112  << std::setfill('0')
113  << std::right
114  << a_step;
115 
117  (file_name_gmv.str(), es);
118  }
119 #endif
120 
121 #ifdef LIBMESH_HAVE_EXODUS_API
122  if (param.output_exodus)
123  {
124  MeshBase & mesh = es.get_mesh();
125 
126  // We write out one file per adaptive step. The files are named in
127  // the following way:
128  // foo.e
129  // foo.e-s002
130  // foo.e-s003
131  // ...
132  // so that, if you open the first one with Paraview, it actually
133  // opens the entire sequence of adapted files.
134  std::ostringstream file_name_exodus;
135 
136  file_name_exodus << solution_type << ".e";
137  if (a_step > 0)
138  file_name_exodus << "-s"
139  << std::setw(3)
140  << std::setfill('0')
141  << std::right
142  << a_step + 1;
143 
144  // We write each adaptive step as a pseudo "time" step, where the
145  // time simply matches the (1-based) adaptive step we are on.
146  ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
147  es,
148  1,
149  /*time=*/a_step + 1);
150  }
151 #endif
152 }
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
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
MeshBase & mesh
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
This is the MeshBase class.
Definition: mesh_base.h:68
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 &)
const MeshBase & get_mesh() const