libMesh
Functions
adjoints_ex1.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< ErrorEstimatorbuild_error_estimator (FEMParameters &param, QoISet &qois)
 
int main (int argc, char **argv)
 

Function Documentation

UniquePtr<ErrorEstimator> build_error_estimator ( FEMParameters param,
QoISet qois 
)

Definition at line 239 of file adjoints_ex1.C.

References libMesh::AdjointResidualErrorEstimator::dual_error_estimator(), libMesh::AdjointResidualErrorEstimator::error_plot_suffix, libMesh::H1_SEMINORM, FEMParameters::indicator_type, libMesh::out, FEMParameters::patch_reuse, libMesh::AdjointResidualErrorEstimator::primal_error_estimator(), libMesh::AdjointResidualErrorEstimator::qoi_set(), and libMesh::PatchRecoveryErrorEstimator::set_patch_reuse().

Referenced by main().

241 {
242  if (param.indicator_type == "kelly")
243  {
244  libMesh::out << "Using Kelly Error Estimator" << std::endl;
245 
247  }
248  else if (param.indicator_type == "adjoint_residual")
249  {
250  libMesh::out << "Using Adjoint Residual Error Estimator with Patch Recovery Weights" << std::endl;
251 
252  AdjointResidualErrorEstimator * adjoint_residual_estimator = new AdjointResidualErrorEstimator;
253 
254  adjoint_residual_estimator->qoi_set() = qois;
255 
256  adjoint_residual_estimator->error_plot_suffix = "error.gmv";
257 
259  adjoint_residual_estimator->primal_error_estimator().reset(p1);
260 
262  adjoint_residual_estimator->dual_error_estimator().reset(p2);
263 
264  adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
265  p1->set_patch_reuse(param.patch_reuse);
266 
267  adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
268  p2->set_patch_reuse(param.patch_reuse);
269 
270  return UniquePtr<ErrorEstimator>(adjoint_residual_estimator);
271  }
272  else
273  libmesh_error_msg("Unknown indicator_type = " << param.indicator_type);
274 }
This class implements a goal oriented error indicator, by weighting residual-based estimates from the...
QoISet & qoi_set()
Access to the QoISet (default: weight all QoIs equally) to use when computing errors.
std::string indicator_type
UniquePtr< ErrorEstimator > & primal_error_estimator()
Access to the "subestimator" (default: PatchRecovery) to use on the primal/forward solution...
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
std::string error_plot_suffix
To aid in investigating error estimator behavior, set this string to a suffix with which to plot (pre...
This class implements the Kelly error indicator which is based on the flux jumps between elements...
This class implements the Patch Recovery error indicator.
OStreamProxy out
UniquePtr< ErrorEstimator > & dual_error_estimator()
Access to the "subestimator" (default: PatchRecovery) to use on the dual/adjoint solution.
UniquePtr<MeshRefinement> build_mesh_refinement ( MeshBase mesh,
FEMParameters param 
)

Definition at line 216 of file adjoints_ex1.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().

218 {
219  MeshRefinement * mesh_refinement = new MeshRefinement(mesh);
220  mesh_refinement->coarsen_by_parents() = true;
221  mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
222  mesh_refinement->nelem_target() = param.nelem_target;
223  mesh_refinement->refine_fraction() = param.refine_fraction;
224  mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
225  mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
226 
227  return UniquePtr<MeshRefinement>(mesh_refinement);
228 }
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 277 of file adjoints_ex1.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_error_estimator(), build_mesh_refinement(), libMesh::LibMeshInit::comm(), libMesh::default_solver_package(), 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().

278 {
279  // Initialize libMesh.
280  LibMeshInit init (argc, argv);
281 
282  // This example requires a linear solver package.
283  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
284  "--enable-petsc, --enable-trilinos, or --enable-eigen");
285 
286  // Skip adaptive examples on a non-adaptive libMesh build
287 #ifndef LIBMESH_ENABLE_AMR
288  libmesh_example_requires(false, "--enable-amr");
289 #else
290 
291  libMesh::out << "Started " << argv[0] << std::endl;
292 
293  // Make sure the general input file exists, and parse it
294  {
295  std::ifstream i("general.in");
296  if (!i)
297  libmesh_error_msg('[' << init.comm().rank() << "] Can't find general.in; exiting early.");
298  }
299  GetPot infile("general.in");
300 
301  // Read in parameters from the input file
302  FEMParameters param(init.comm());
303  param.read(infile);
304 
305  // Skip this default-2D example if libMesh was compiled as 1D-only.
306  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
307 
308  // Create a mesh, with dimension to be overridden later, distributed
309  // across the default MPI communicator.
310  Mesh mesh(init.comm());
311 
312  // And an object to refine it
313  UniquePtr<MeshRefinement> mesh_refinement =
314  build_mesh_refinement(mesh, param);
315 
316  // And an EquationSystems to run on it
317  EquationSystems equation_systems (mesh);
318 
319  libMesh::out << "Reading in and building the mesh" << std::endl;
320 
321  // Read in the mesh
322  mesh.read(param.domainfile.c_str());
323  // Make all the elements of the mesh second order so we can compute
324  // with a higher order basis
326 
327  // Create a mesh refinement object to do the initial uniform refinements
328  // on the coarse grid read in from lshaped.xda
329  MeshRefinement initial_uniform_refinements(mesh);
330  initial_uniform_refinements.uniformly_refine(param.coarserefinements);
331 
332  libMesh::out << "Building system" << std::endl;
333 
334  // Build the FEMSystem
335  LaplaceSystem & system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
336 
337  // Set its parameters
338  set_system_parameters(system, param);
339 
340  libMesh::out << "Initializing systems" << std::endl;
341 
342  equation_systems.init ();
343 
344  // Print information about the mesh and system to the screen.
345  mesh.print_info();
346  equation_systems.print_info();
347  LinearSolver<Number> * linear_solver = system.get_linear_solver();
348 
349  {
350  // Adaptively solve the timestep
351  unsigned int a_step = 0;
352  for (; a_step != param.max_adaptivesteps; ++a_step)
353  {
354  // We can't adapt to both a tolerance and a
355  // target mesh size
356  if (param.global_tolerance != 0.)
357  libmesh_assert_equal_to (param.nelem_target, 0);
358  // If we aren't adapting to a tolerance we need a
359  // target mesh size
360  else
361  libmesh_assert_greater (param.nelem_target, 0);
362 
363  linear_solver->reuse_preconditioner(false);
364 
365  // Solve the forward problem
366  system.solve();
367 
368  // Write out the computed primal solution
369  write_output(equation_systems, a_step, "primal", param);
370 
371  // Get a pointer to the primal solution vector
372  NumericVector<Number> & primal_solution = *system.solution;
373 
374  // Declare a QoISet object, we need this object to set weights for our QoI error contributions
375  QoISet qois;
376 
377  // Declare a qoi_indices vector, each index will correspond to a QoI
378  std::vector<unsigned int> qoi_indices;
379  qoi_indices.push_back(0);
380  qoi_indices.push_back(1);
381  qois.add_indices(qoi_indices);
382 
383  // Set weights for each index, these will weight the contribution of each QoI in the final error
384  // estimate to be used for flagging elements for refinement
385  qois.set_weight(0, 0.5);
386  qois.set_weight(1, 0.5);
387 
388  // Make sure we get the contributions to the adjoint RHS from the sides
389  system.assemble_qoi_sides = true;
390 
391  // We are about to solve the adjoint system, but before we do this we see the same preconditioner
392  // flag to reuse the preconditioner from the forward solver
393  linear_solver->reuse_preconditioner(param.reuse_preconditioner);
394 
395  // Solve the adjoint system. This takes the transpose of the stiffness matrix and then
396  // solves the resulting system
397  system.adjoint_solve();
398 
399  // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
400  system.set_adjoint_already_solved(true);
401 
402  // Get a pointer to the solution vector of the adjoint problem for QoI 0
403  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
404 
405  // Swap the primal and dual solutions so we can write out the adjoint solution
406  primal_solution.swap(dual_solution_0);
407  write_output(equation_systems, a_step, "adjoint_0", param);
408 
409  // Swap back
410  primal_solution.swap(dual_solution_0);
411 
412  // Get a pointer to the solution vector of the adjoint problem for QoI 0
413  NumericVector<Number> & dual_solution_1 = system.get_adjoint_solution(1);
414 
415  // Swap again
416  primal_solution.swap(dual_solution_1);
417  write_output(equation_systems, a_step, "adjoint_1", param);
418 
419  // Swap back again
420  primal_solution.swap(dual_solution_1);
421 
422  libMesh::out << "Adaptive step "
423  << a_step
424  << ", we have "
425  << mesh.n_active_elem()
426  << " active elements and "
427  << equation_systems.n_active_dofs()
428  << " active dofs."
429  << std::endl;
430 
431  // Postprocess, compute the approximate QoIs and write them out to the console
432  libMesh::out << "Postprocessing: " << std::endl;
433  system.postprocess_sides = true;
434  system.postprocess();
435  Number QoI_0_computed = system.get_QoI_value("computed", 0);
436  Number QoI_0_exact = system.get_QoI_value("exact", 0);
437  Number QoI_1_computed = system.get_QoI_value("computed", 1);
438  Number QoI_1_exact = system.get_QoI_value("exact", 1);
439 
440  libMesh::out << "The relative error in QoI 0 is "
441  << std::setprecision(17)
442  << std::abs(QoI_0_computed - QoI_0_exact) / std::abs(QoI_0_exact)
443  << std::endl;
444 
445  libMesh::out << "The relative error in QoI 1 is "
446  << std::setprecision(17)
447  << std::abs(QoI_1_computed - QoI_1_exact) / std::abs(QoI_1_exact)
448  << std::endl
449  << std::endl;
450 
451  // Now we construct the data structures for the mesh refinement process
452  ErrorVector error;
453 
454  // Build an error estimator object
455  UniquePtr<ErrorEstimator> error_estimator =
456  build_error_estimator(param, qois);
457 
458  // Estimate the error in each element using the Adjoint Residual or Kelly error estimator
459  error_estimator->estimate_error(system, error);
460 
461  // We have to refine either based on reaching an error tolerance or
462  // a number of elements target, which should be verified above
463  // Otherwise we flag elements by error tolerance or nelem target
464 
465  // Uniform refinement
466  if (param.refine_uniformly)
467  {
468  mesh_refinement->uniformly_refine(1);
469  }
470  // Adaptively refine based on reaching an error tolerance
471  else if (param.global_tolerance >= 0. && param.nelem_target == 0.)
472  {
473  mesh_refinement->flag_elements_by_error_tolerance (error);
474 
475  mesh_refinement->refine_and_coarsen_elements();
476  }
477  // Adaptively refine based on reaching a target number of elements
478  else
479  {
480  if (mesh.n_active_elem() >= param.nelem_target)
481  {
482  libMesh::out << "We reached the target number of elements." << std::endl << std::endl;
483  break;
484  }
485 
486  mesh_refinement->flag_elements_by_nelem_target (error);
487 
488  mesh_refinement->refine_and_coarsen_elements();
489  }
490 
491  // Dont forget to reinit the system after each adaptive refinement !
492  equation_systems.reinit();
493 
494  libMesh::out << "Refined mesh to "
495  << mesh.n_active_elem()
496  << " active elements and "
497  << equation_systems.n_active_dofs()
498  << " active dofs."
499  << std::endl;
500  }
501 
502  // Do one last solve if necessary
503  if (a_step == param.max_adaptivesteps)
504  {
505  linear_solver->reuse_preconditioner(false);
506  system.solve();
507 
508  write_output(equation_systems, a_step, "primal", param);
509 
510  NumericVector<Number> & primal_solution = *system.solution;
511 
512  QoISet qois;
513  std::vector<unsigned int> qoi_indices;
514 
515  qoi_indices.push_back(0);
516  qoi_indices.push_back(1);
517  qois.add_indices(qoi_indices);
518 
519  qois.set_weight(0, 0.5);
520  qois.set_weight(1, 0.5);
521 
522  system.assemble_qoi_sides = true;
523  linear_solver->reuse_preconditioner(param.reuse_preconditioner);
524  system.adjoint_solve();
525 
526  // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
527  system.set_adjoint_already_solved(true);
528 
529  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
530 
531  primal_solution.swap(dual_solution_0);
532  write_output(equation_systems, a_step, "adjoint_0", param);
533 
534  primal_solution.swap(dual_solution_0);
535 
536  NumericVector<Number> & dual_solution_1 = system.get_adjoint_solution(1);
537 
538  primal_solution.swap(dual_solution_1);
539  write_output(equation_systems, a_step, "adjoint_1", param);
540 
541  primal_solution.swap(dual_solution_1);
542 
543  libMesh::out << "Adaptive step "
544  << a_step
545  << ", we have "
546  << mesh.n_active_elem()
547  << " active elements and "
548  << equation_systems.n_active_dofs()
549  << " active dofs."
550  << std::endl;
551 
552  libMesh::out << "Postprocessing: " << std::endl;
553  system.postprocess_sides = true;
554  system.postprocess();
555 
556  Number QoI_0_computed = system.get_QoI_value("computed", 0);
557  Number QoI_0_exact = system.get_QoI_value("exact", 0);
558  Number QoI_1_computed = system.get_QoI_value("computed", 1);
559  Number QoI_1_exact = system.get_QoI_value("exact", 1);
560 
561  libMesh::out << "The relative error in QoI 0 is "
562  << std::setprecision(17)
563  << std::abs(QoI_0_computed - QoI_0_exact) / std::abs(QoI_0_exact)
564  << std::endl;
565 
566  libMesh::out << "The relative error in QoI 1 is "
567  << std::setprecision(17)
568  << std::abs(QoI_1_computed - QoI_1_exact) / std::abs(QoI_1_exact)
569  << std::endl
570  << std::endl;
571 
572  // Hard coded asserts to ensure that the actual numbers we are getting are what they should be
573  libmesh_assert_less(std::abs(QoI_0_computed - QoI_0_exact)/std::abs(QoI_0_exact), 4.e-5);
574  libmesh_assert_less(std::abs(QoI_1_computed - QoI_1_exact)/std::abs(QoI_1_exact), 1.e-4);
575  }
576  }
577 
578  libMesh::err << '[' << mesh.processor_id()
579  << "] Completing output."
580  << std::endl;
581 
582 #endif // #ifndef LIBMESH_ENABLE_AMR
583 
584  // All done.
585  return 0;
586 }
OStreamProxy err
double abs(double a)
This is the EquationSystems class.
virtual dof_id_type n_active_elem() const =0
void write_output(EquationSystems &es, unsigned int a_step, std::string solution_type, FEMParameters &param)
Definition: adjoints_ex1.C:97
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
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
void set_system_parameters(LaplaceSystem &system, FEMParameters &param)
Definition: adjoints_ex1.C:161
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.
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 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
UniquePtr< ErrorEstimator > build_error_estimator(FEMParameters &param, QoISet &qois)
Definition: adjoints_ex1.C:239
UniquePtr< MeshRefinement > build_mesh_refinement(MeshBase &mesh, FEMParameters &param)
Definition: adjoints_ex1.C:216
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.
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 161 of file adjoints_ex1.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, FEMParameters::solver_verbose, libMesh::DifferentiableSystem::time_solver, libMesh::DiffSolver::verbose, FEMParameters::verify_analytic_jacobians, and libMesh::FEMSystem::verify_analytic_jacobians.

Referenced by main().

163 {
164  // Use analytical jacobians?
165  system.analytic_jacobians() = param.analytic_jacobians;
166 
167  // Verify analytic jacobians against numerical ones?
169 
170  // Use the prescribed FE type
171  system.fe_family() = param.fe_family[0];
172  system.fe_order() = param.fe_order[0];
173 
174  // More desperate debugging options
176  system.print_solutions = param.print_solutions;
178  system.print_residuals = param.print_residuals;
180  system.print_jacobians = param.print_jacobians;
181 
182  // No transient time solver
183  system.time_solver =
184  UniquePtr<TimeSolver>(new SteadySolver(system));
185 
186  // Nonlinear solver options
187  {
188  NewtonSolver * solver = new NewtonSolver(system);
189  system.time_solver->diff_solver() = UniquePtr<DiffSolver>(solver);
190 
191  solver->quiet = param.solver_quiet;
192  solver->verbose = param.solver_verbose;
194  solver->minsteplength = param.min_step_length;
199  if (system.time_solver->reduce_deltat_on_diffsolver_failure)
200  {
201  solver->continue_after_max_iterations = true;
202  solver->continue_after_backtrack_failure = true;
203  }
204 
205  // And the linear solver options
209  }
210 }
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 verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
Definition: diff_solver.h:168
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 97 of file adjoints_ex1.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().

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