libMesh
Functions
adjoints_ex5.C File Reference

Go to the source code of this file.

Functions

void write_output (EquationSystems &es, unsigned int t_step, std::string solution_type, FEMParameters &param)
 
void set_system_parameters (HeatSystem &system, FEMParameters &param)
 
int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 290 of file adjoints_ex5.C.

References std::abs(), libMesh::System::add_vector(), adjoint_initial_grad(), adjoint_initial_value(), libMesh::DifferentiableSystem::adjoint_solve(), libMesh::MeshTools::Generation::build_square(), libMesh::System::calculate_norm(), libMesh::default_solver_package(), libMesh::DifferentiableSystem::deltat, libMesh::err, finish_initialization(), libMesh::System::get_adjoint_solution(), libMesh::System::get_vector(), libMesh::GHOSTED, libMesh::H1, libMesh::TriangleWrapper::init(), initial_grad(), initial_value(), libMesh::libmesh_ignore(), mesh, libMesh::out, libMesh::PETSC_SOLVERS, libMesh::System::project_solution(), libMesh::System::project_vector(), libMesh::QUAD4, FEMParameters::read(), read_initial_parameters(), libMesh::System::set_adjoint_already_solved(), set_system_parameters(), libMesh::System::set_vector_preservation(), libMesh::System::solution, libMesh::FEMSystem::solve(), libMesh::NumericVector< T >::swap(), libMesh::System::time, libMesh::DifferentiableSystem::time_solver, libMesh::TRI3, and write_output().

291 {
292  // Skip adaptive examples on a non-adaptive libMesh build
293 #ifndef LIBMESH_ENABLE_AMR
294  libmesh_ignore(argc);
295  libmesh_ignore(argv);
296  libmesh_example_requires(false, "--enable-amr");
297 #else
298  // Skip this 2D example if libMesh was compiled as 1D-only.
299  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
300 
301  // Initialize libMesh.
302  LibMeshInit init (argc, argv);
303 
304  // This doesn't converge with Trilinos for some reason...
305  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
306 
307  libMesh::out << "Started " << argv[0] << std::endl;
308 
309  // Make sure the general input file exists, and parse it
310  {
311  std::ifstream i("general.in");
312  if (!i)
313  libmesh_error_msg('[' << init.comm().rank() << "] Can't find general.in; exiting early.");
314  }
315  GetPot infile("general.in");
316 
317  // Read in parameters from the input file
318  FEMParameters param(init.comm());
319  param.read(infile);
320 
321  // Create a mesh with the given dimension, distributed
322  // across the default MPI communicator.
323  Mesh mesh(init.comm(), param.dimension);
324 
325  // And an object to refine it
326  UniquePtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
327 
328  // And an EquationSystems to run on it
329  EquationSystems equation_systems (mesh);
330 
331  libMesh::out << "Building mesh" << std::endl;
332 
333  // Build a unit square
334  ElemType elemtype;
335 
336  if (param.elementtype == "tri" ||
337  param.elementtype == "unstructured")
338  elemtype = TRI3;
339  else
340  elemtype = QUAD4;
341 
342  MeshTools::Generation::build_square (mesh, param.coarsegridx, param.coarsegridy,
343  param.domain_xmin, param.domain_xmin + param.domain_edge_width,
344  param.domain_ymin, param.domain_ymin + param.domain_edge_length,
345  elemtype);
346 
347  libMesh::out << "Building system" << std::endl;
348 
349  HeatSystem & system = equation_systems.add_system<HeatSystem> ("HeatSystem");
350 
351  set_system_parameters(system, param);
352 
353  libMesh::out << "Initializing systems" << std::endl;
354 
355  // Initialize the system
356  equation_systems.init ();
357 
358  // Refine the grid again if requested
359  for (unsigned int i=0; i != param.extrarefinements; ++i)
360  {
361  mesh_refinement->uniformly_refine(1);
362  equation_systems.reinit();
363  }
364 
365  libMesh::out << "Setting primal initial conditions" << std::endl;
366 
368 
370  equation_systems.parameters);
371 
372  // Output the H1 norm of the initial conditions
373  libMesh::out << "|U("
374  << system.time
375  << ")|= "
376  << system.calculate_norm(*system.solution, 0, H1)
377  << std::endl
378  << std::endl;
379 
380  // Add an adjoint vector, this will be computed after the forward
381  // time stepping is complete
382  //
383  // Tell the library not to save adjoint solutions during the forward
384  // solve
385  //
386  // Tell the library not to project this vector, and hence, memory
387  // solution history to not save it.
388  //
389  // Make this vector ghosted so we can localize it to each element
390  // later.
391  const std::string & adjoint_solution_name = "adjoint_solution0";
392  system.add_vector("adjoint_solution0", false, GHOSTED);
393 
394  // Close up any resources initial.C needed
396 
397  // Plot the initial conditions
398  write_output(equation_systems, 0, "primal", param);
399 
400  // Print information about the mesh and system to the screen.
401  mesh.print_info();
402  equation_systems.print_info();
403 
404  // In optimized mode we catch any solver errors, so that we can
405  // write the proper footers before closing. In debug mode we just
406  // let the exception throw so that gdb can grab it.
407 #ifdef NDEBUG
408  try
409  {
410 #endif
411  // Now we begin the timestep loop to compute the time-accurate
412  // solution of the equations.
413  for (unsigned int t_step=param.initial_timestep;
414  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
415  {
416  // A pretty update message
417  libMesh::out << " Solving time step "
418  << t_step
419  << ", time = "
420  << system.time
421  << std::endl;
422 
423  // Solve the forward problem at time t, to obtain the solution at time t + dt
424  system.solve();
425 
426  // Output the H1 norm of the computed solution
427  libMesh::out << "|U("
428  << system.time + system.deltat
429  << ")|= "
430  << system.calculate_norm(*system.solution, 0, H1)
431  << std::endl;
432 
433  // Advance to the next timestep in a transient problem
434  libMesh::out << "Advancing timestep" << std::endl << std::endl;
435  system.time_solver->advance_timestep();
436 
437  // Write out this timestep
438  write_output(equation_systems, t_step+1, "primal", param);
439  }
440  // End timestep loop
441 
443 
444  // Now we will solve the backwards in time adjoint problem
445  libMesh::out << std::endl << "Solving the adjoint problem" << std::endl;
446 
447  // We need to tell the library that it needs to project the adjoint, so
448  // MemorySolutionHistory knows it has to save it
449 
450  // Tell the library to project the adjoint vector, and hence, memory solution history to
451  // save it
452  system.set_vector_preservation(adjoint_solution_name, true);
453 
454  libMesh::out << "Setting adjoint initial conditions Z("
455  << system.time
456  << ")"
457  <<std::endl;
458 
459  // Need to call adjoint_advance_timestep once for the initial condition setup
460  libMesh::out<<"Retrieving solutions at time t="<<system.time<<std::endl;
461  system.time_solver->adjoint_advance_timestep();
462 
463  // Output the H1 norm of the retrieved solutions (u^i and u^i+1)
464  libMesh::out << "|U("
465  << system.time + system.deltat
466  << ")|= "
467  << system.calculate_norm(*system.solution, 0, H1)
468  << std::endl;
469 
470  libMesh::out << "|U("
471  << system.time
472  << ")|= "
473  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
474  << std::endl;
475 
476  // The first thing we have to do is to apply the adjoint initial
477  // condition. The user should supply these. Here they are specified
478  // in the functions adjoint_initial_value and adjoint_initial_gradient
481  equation_systems.parameters,
482  system.get_adjoint_solution(0));
483 
484  // Since we have specified an adjoint solution for the current
485  // time (T), set the adjoint_already_solved boolean to true, so
486  // we dont solve unnecessarily in the adjoint sensitivity method
487  system.set_adjoint_already_solved(true);
488 
489  libMesh::out << "|Z("
490  << system.time
491  << ")|= "
492  << system.calculate_norm(system.get_adjoint_solution(), 0, H1)
493  << std::endl
494  << std::endl;
495 
496  write_output(equation_systems, param.n_timesteps, "dual", param);
497 
498  // Now that the adjoint initial condition is set, we will start the
499  // backwards in time adjoint integration
500 
501  // For loop stepping backwards in time
502  for (unsigned int t_step=param.initial_timestep;
503  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
504  {
505  //A pretty update message
506  libMesh::out << " Solving adjoint time step "
507  << t_step
508  << ", time = "
509  << system.time
510  << std::endl;
511 
512  // The adjoint_advance_timestep function calls the retrieve
513  // function of the memory_solution_history class via the
514  // memory_solution_history object we declared earlier. The
515  // retrieve function sets the system primal vectors to their
516  // values at the current timestep.
517  libMesh::out << "Retrieving solutions at time t=" << system.time << std::endl;
518  system.time_solver->adjoint_advance_timestep();
519 
520  // Output the H1 norm of the retrieved solution
521  libMesh::out << "|U("
522  << system.time + system.deltat
523  << ")|= "
524  << system.calculate_norm(*system.solution, 0, H1)
525  << std::endl;
526 
527  libMesh::out << "|U("
528  << system.time
529  << ")|= "
530  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
531  << std::endl;
532 
533  system.set_adjoint_already_solved(false);
534 
535  system.adjoint_solve();
536 
537  // Now that we have solved the adjoint, set the
538  // adjoint_already_solved boolean to true, so we dont solve
539  // unnecessarily in the error estimator
540  system.set_adjoint_already_solved(true);
541 
542  libMesh::out << "|Z("
543  << system.time
544  << ")|= "
545  << system.calculate_norm(system.get_adjoint_solution(), 0, H1)
546  << std::endl
547  << std::endl;
548 
549  // Get a pointer to the primal solution vector
550  NumericVector<Number> & primal_solution = *system.solution;
551 
552  // Get a pointer to the solution vector of the adjoint problem for QoI 0
553  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
554 
555  // Swap the primal and dual solutions so we can write out the adjoint solution
556  primal_solution.swap(dual_solution_0);
557 
558  write_output(equation_systems, param.n_timesteps - (t_step + 1), "dual", param);
559 
560  // Swap back
561  primal_solution.swap(dual_solution_0);
562  }
563  // End adjoint timestep loop
564 
565  // Now that we have computed both the primal and adjoint solutions, we compute the sensitivities to the parameter p
566  // dQ/dp = partialQ/partialp - partialR/partialp
567  // partialQ/partialp = (Q(p+dp) - Q(p-dp))/(2*dp), this is not supported by the library yet
568  // partialR/partialp = (R(u,z;p+dp) - R(u,z;p-dp))/(2*dp), where
569  // R(u,z;p+dp) = int_{0}^{T} f(z;p+dp) - <partialu/partialt, z>(p+dp) - <g(u),z>(p+dp)
570  // To do this we need to step forward in time, and compute the perturbed R at each time step and accumulate it
571  // Then once all time steps are over, we can compute (R(u,z;p+dp) - R(u,z;p-dp))/(2*dp)
572 
573  // Now we begin the timestep loop to compute the time-accurate
574  // adjoint sensitivities
575  for (unsigned int t_step=param.initial_timestep;
576  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
577  {
578  // A pretty update message
579  libMesh::out << "Retrieving "
580  << t_step
581  << ", time = "
582  << system.time
583  << std::endl;
584 
585  // Retrieve the primal and adjoint solutions at the current timestep
586  system.time_solver->retrieve_timestep();
587 
588  libMesh::out << "|U("
589  << system.time + system.deltat
590  << ")|= "
591  << system.calculate_norm(*system.solution, 0, H1)
592  << std::endl;
593 
594  libMesh::out << "|U("
595  << system.time
596  << ")|= "
597  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
598  << std::endl;
599 
600  libMesh::out << "|Z("
601  << system.time
602  << ")|= "
603  << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
604  << std::endl
605  << std::endl;
606 
607  // Call the postprocess function which we have overloaded to compute
608  // accumulate the perturbed residuals
609  dynamic_cast<HeatSystem &>(system).perturb_accumulate_residuals(dynamic_cast<HeatSystem &>(system).get_parameter_vector());
610 
611  // Move the system time forward (retrieve_timestep does not do this)
612  system.time += system.deltat;
613  }
614 
615  // A pretty update message
616  libMesh::out << "Retrieving final time = "
617  << system.time
618  << std::endl;
619 
620  // Retrieve the primal and adjoint solutions at the current timestep
621  system.time_solver->retrieve_timestep();
622 
623  libMesh::out << "|U("
624  << system.time + system.deltat
625  << ")|= "
626  << system.calculate_norm(*system.solution, 0, H1)
627  << std::endl;
628 
629  libMesh::out << "|U("
630  << system.time
631  << ")|= "
632  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
633  << std::endl;
634 
635  libMesh::out << "|Z("
636  << system.time
637  << ")|= "
638  << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
639  << std::endl
640  << std::endl;
641 
642  // Call the postprocess function which we have overloaded to compute
643  // accumulate the perturbed residuals
644  dynamic_cast<HeatSystem &>(system).perturb_accumulate_residuals(dynamic_cast<HeatSystem &>(system).get_parameter_vector());
645 
646  // Now that we computed the accumulated, perturbed residuals, we can compute the
647  // approximate sensitivity
648  Number sensitivity_0_0 = (dynamic_cast<HeatSystem &>(system)).compute_final_sensitivity();
649 
650  // Print it out
651  libMesh::out << "Sensitivity of QoI 0 w.r.t parameter 0 is: "
652  << sensitivity_0_0
653  << std::endl;
654 
655  // Hard coded assert to ensure that the actual numbers we are
656  // getting are what they should be
657  // The 2e-4 tolerance is chosen to ensure success even with
658  // 32-bit floats
659  libmesh_assert_less(std::abs(sensitivity_0_0 - (-5.37173)), 2.e-4);
660 
661 #ifdef NDEBUG
662  }
663  catch (...)
664  {
665  libMesh::err << '[' << mesh.processor_id()
666  << "] Caught exception; exiting early." << std::endl;
667  }
668 #endif
669 
670  libMesh::err << '[' << mesh.processor_id()
671  << "] Completing output."
672  << std::endl;
673 
674  // All done.
675  return 0;
676 
677 #endif // LIBMESH_ENABLE_AMR
678 }
OStreamProxy err
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1545
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
void write_output(EquationSystems &es, unsigned int t_step, std::string solution_type, FEMParameters &param)
Definition: adjoints_ex5.C:111
void finish_initialization()
Definition: initial.C:11
void set_adjoint_already_solved(bool setting)
Setter for the adjoint_already_solved boolean.
Definition: system.h:378
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
ElemType
Defines an enum for geometric element types.
MeshBase & mesh
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.
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr) const
Projects arbitrary functions onto the current solution.
Number adjoint_initial_value(const Point &p, const Parameters &, const std::string &, const std::string &)
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
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
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
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
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
Real calculate_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type, std::set< unsigned int > *skip_dimensions=libmesh_nullptr) const
Definition: system.C:1383
void libmesh_ignore(const T &)
void set_system_parameters(HeatSystem &system, FEMParameters &param)
Definition: adjoints_ex5.C:171
void set_vector_preservation(const std::string &vec_name, bool preserve)
Allows one to set the boolean controlling whether the vector identified by vec_name should be "preser...
Definition: system.C:887
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
OStreamProxy out
Gradient initial_grad(const Point &, const Parameters &, const std::string &, const std::string &)
Definition: initial.C:30
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
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
Projects arbitrary functions onto a vector of degree of freedom values for the current system...
void read_initial_parameters()
Definition: initial.C:7
Number initial_value(const Point &, const Parameters &, const std::string &, const std::string &)
Definition: initial.C:18
Gradient adjoint_initial_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
void set_system_parameters ( HeatSystem system,
FEMParameters param 
)

Definition at line 171 of file adjoints_ex5.C.

References libMesh::DofMap::add_dirichlet_boundary(), HeatSystem::analytic_jacobians(), FEMParameters::analytic_jacobians, FEMParameters::deltat, libMesh::DifferentiableSystem::deltat, FEMParameters::deltat_reductions, FEMParameters::dirichlet_condition_variables, FEMParameters::dirichlet_conditions, FEMParameters::extra_quadrature_order, libMesh::System::extra_quadrature_order, HeatSystem::fe_family(), FEMParameters::fe_family, HeatSystem::fe_order(), FEMParameters::fe_order, libMesh::System::get_dof_map(), FEMParameters::initial_linear_tolerance, FEMParameters::linear_tolerance_multiplier, FEMParameters::max_linear_iterations, FEMParameters::max_nonlinear_iterations, FEMParameters::min_step_length, FEMParameters::minimum_linear_tolerance, FEMParameters::numerical_jacobian_h, libMesh::FEMSystem::numerical_jacobian_h, libMesh::out, FEMParameters::print_element_jacobians, libMesh::DifferentiableSystem::print_element_jacobians, FEMParameters::print_element_residuals, libMesh::DifferentiableSystem::print_element_residuals, 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, FEMParameters::relative_residual_tolerance, FEMParameters::relative_step_tolerance, FEMParameters::require_residual_reduction, libMesh::solver, FEMParameters::solver_quiet, FEMParameters::solver_verbose, libMesh::DifferentiableSystem::time_solver, FEMParameters::time_solver_quiet, FEMParameters::timesolver_core, FEMParameters::timesolver_theta, FEMParameters::transient, FEMParameters::use_petsc_snes, FEMParameters::verify_analytic_jacobians, and libMesh::FEMSystem::verify_analytic_jacobians.

Referenced by main().

173 {
174  // Use the prescribed FE type
175  system.fe_family() = param.fe_family[0];
176  system.fe_order() = param.fe_order[0];
177 
178  // Use analytical jacobians?
179  system.analytic_jacobians() = param.analytic_jacobians;
180 
181  // Verify analytic jacobians against numerical ones?
184 
185  // More desperate debugging options
187  system.print_solutions = param.print_solutions;
189  system.print_residuals = param.print_residuals;
191  system.print_jacobians = param.print_jacobians;
194 
195  // Solve this as a time-dependent or steady system
196  if (param.transient)
197  {
198  UnsteadySolver *innersolver;
199  if (param.timesolver_core == "euler")
200  {
201  EulerSolver *eulersolver =
202  new EulerSolver(system);
203 
204  eulersolver->theta = param.timesolver_theta;
205  innersolver = eulersolver;
206  }
207  else
208  libmesh_error_msg("This example (and unsteady adjoints in libMesh) only support Backward Euler and explicit methods.");
209 
210  system.time_solver =
211  UniquePtr<TimeSolver>(innersolver);
212  }
213  else
214  system.time_solver =
215  UniquePtr<TimeSolver>(new SteadySolver(system));
216 
217  // The Memory Solution History object we will set the system SolutionHistory object to
218  MemorySolutionHistory heatsystem_solution_history(system);
219  system.time_solver->set_solution_history(heatsystem_solution_history);
220 
221  system.time_solver->reduce_deltat_on_diffsolver_failure =
222  param.deltat_reductions;
223  system.time_solver->quiet = param.time_solver_quiet;
224 
225  // Create any Dirichlet boundary conditions
226  typedef
227  std::map<boundary_id_type, FunctionBase<Number> *>::
228  const_iterator Iter;
229 
230  for (Iter i = param.dirichlet_conditions.begin();
231  i != param.dirichlet_conditions.end(); ++i)
232  {
233  boundary_id_type b = i->first;
234  FunctionBase<Number> *f = i->second;
235  std::set<boundary_id_type> bdys; bdys.insert(b);
236 
237  system.get_dof_map().add_dirichlet_boundary(DirichletBoundary(bdys,
239  f));
240 
241  libMesh::out << "Added Dirichlet boundary " << b << " for variables ";
242  for (std::size_t vi=0; vi != param.dirichlet_condition_variables[b].size(); ++vi)
244  libMesh::out << std::endl;
245  }
246 
247  // Set the time stepping options
248  system.deltat = param.deltat;
249 
250  // And the integration options
252 
253  // And the nonlinear solver options
254  if (param.use_petsc_snes)
255  {
256 #ifdef LIBMESH_HAVE_PETSC
257  PetscDiffSolver *solver = new PetscDiffSolver(system);
258  system.time_solver->diff_solver() = UniquePtr<DiffSolver>(solver);
259 #else
260  libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
261 #endif
262  }
263  else
264  {
265  NewtonSolver *solver = new NewtonSolver(system);
266  system.time_solver->diff_solver() = UniquePtr<DiffSolver>(solver);
267 
268  solver->quiet = param.solver_quiet;
269  solver->verbose = param.solver_verbose;
270  solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
271  solver->minsteplength = param.min_step_length;
272  solver->relative_step_tolerance = param.relative_step_tolerance;
273  solver->relative_residual_tolerance = param.relative_residual_tolerance;
274  solver->require_residual_reduction = param.require_residual_reduction;
275  solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
276  if (system.time_solver->reduce_deltat_on_diffsolver_failure)
277  {
278  solver->continue_after_max_iterations = true;
279  solver->continue_after_backtrack_failure = true;
280  }
281 
282  // And the linear solver options
283  solver->max_linear_iterations = param.max_linear_iterations;
284  solver->initial_linear_tolerance = param.initial_linear_tolerance;
285  solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
286  }
287 }
bool print_solution_norms
libMesh::Real deltat
Definition: femparameters.h:38
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
bool print_element_residuals
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 print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Definition: diff_system.h:329
int extra_quadrature_order
A member int that can be employed to indicate increased or reduced quadrature order.
Definition: system.h:1508
Real numerical_jacobian_h
If calculating numeric jacobians is required, the FEMSystem will perturb each solution vector entry b...
Definition: fem_system.h:175
std::string timesolver_core
Definition: femparameters.h:37
bool require_residual_reduction
bool print_residual_norms
libMesh::Real relative_step_tolerance
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:334
std::vector< unsigned int > fe_order
bool print_element_residuals
Set print_element_residuals to true to print each R_elem contribution.
Definition: diff_system.h:344
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
const DofMap & get_dof_map() const
Definition: system.h:2030
unsigned int max_linear_iterations
int8_t boundary_id_type
Definition: id_types.h:51
bool print_residual_norms
Set print_residual_norms to true to print |F| whenever it is assembled.
Definition: diff_system.h:319
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
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
int extra_quadrature_order
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
bool print_jacobian_norms
OStreamProxy out
libMesh::Real min_step_length
std::string & fe_family()
Definition: heatsystem.h:46
libMesh::Real timesolver_theta
Definition: femparameters.h:38
bool print_element_jacobians
Set print_element_jacobians to true to print each J_elem contribution.
Definition: diff_system.h:349
bool print_element_jacobians
unsigned int max_nonlinear_iterations
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
std::map< libMesh::boundary_id_type, std::vector< unsigned int > > dirichlet_condition_variables
Definition: femparameters.h:90
bool time_solver_quiet
bool & analytic_jacobians()
Definition: heatsystem.h:49
libMesh::Real numerical_jacobian_h
std::map< libMesh::boundary_id_type, libMesh::FunctionBase< libMesh::Number > * > dirichlet_conditions
Definition: femparameters.h:87
libMesh::Real linear_tolerance_multiplier
unsigned int & fe_order()
Definition: heatsystem.h:47
unsigned int deltat_reductions
Definition: femparameters.h:36
void write_output ( EquationSystems &  es,
unsigned int  t_step,
std::string  solution_type,
FEMParameters param 
)

Definition at line 111 of file adjoints_ex5.C.

References libMesh::libmesh_ignore(), mesh, FEMParameters::output_exodus, and FEMParameters::output_gmv.

Referenced by main().

115 {
116  // Ignore parameters when there are no output formats available.
117  libmesh_ignore(es);
118  libmesh_ignore(t_step);
119  libmesh_ignore(solution_type);
120  libmesh_ignore(param);
121 
122 #ifdef LIBMESH_HAVE_GMV
123  if (param.output_gmv)
124  {
125  MeshBase & mesh = es.get_mesh();
126 
127  std::ostringstream file_name_gmv;
128  file_name_gmv << solution_type
129  << ".out.gmv."
130  << std::setw(2)
131  << std::setfill('0')
132  << std::right
133  << t_step;
134 
135  GMVIO(mesh).write_equation_systems(file_name_gmv.str(), es);
136  }
137 #endif
138 
139 #ifdef LIBMESH_HAVE_EXODUS_API
140  if (param.output_exodus)
141  {
142  MeshBase & mesh = es.get_mesh();
143 
144  // We write out one file per timestep. The files are named in
145  // the following way:
146  // foo.e
147  // foo.e-s002
148  // foo.e-s003
149  // ...
150  // so that, if you open the first one with Paraview, it actually
151  // opens the entire sequence of adapted files.
152  std::ostringstream file_name_exodus;
153 
154  file_name_exodus << solution_type << ".e";
155  if (t_step > 0)
156  file_name_exodus << "-s"
157  << std::setw(3)
158  << std::setfill('0')
159  << std::right
160  << t_step + 1;
161 
162  // TODO: Get the current time from the System...
163  ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
164  es,
165  1,
166  /*time=*/t_step + 1);
167  }
168 #endif
169 }
MeshBase & mesh
void libmesh_ignore(const T &)