libMesh
adjoints_ex5.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // Adjoints Example 5 - SolutionHistory (both Memory and File), General Localized Vectors and Unsteady Adjoints.
21 // \author Vikram Garg
22 // \date 2020
23 //
24 // This example showcases the solution storage and retrieval capabilities of libMesh, in the context of
25 // unsteady adjoints. The primary motivation is adjoint sensitivity analysis for unsteady
26 // problems. The PDE we are interested in is the simple 2-d heat
27 // equation:
28 // partial(T)/partial(t) - K Laplacian(T) = 0
29 // with initial condition:
30 // T(x,y;0) = sin(pi*x) sin(pi*y) and boundary conditions:
31 // T(boundary;t) = 0
32 
33 // For these initial and boundary conditions, the exact solution
34 // u = exp(-K 2*pi^2 t) * sin(pi*x) * sin(pi*y)
35 
36 // We specify our Quantity of Interest (QoI) as
37 // Q(u) = int_{domain} u(x,y;1) sin(pi*x) sin(pi*y) dx dy, and
38 // are interested in computing the sensitivity dQ/dK
39 
40 // The exact value of this sensitivity is:
41 // dQ/dK = int_{domain} du/dK sin(pi*x) sin(pi*y) dx dy
42 // = int_{domain} (-2*pi^2 * exp(-K pi^2) ) sin^2(pi*x) sin^2(pi*y) dx dy
43 // = (-2*pi^2 * exp(-K 2*pi^2) )/4 = -4.9022 (K = 1.0e-3)
44 
45 // For this QoI, the continuous adjoint problem reads,
46 // -partial(z)/partial(t) - K Laplacian(z) = 0
47 // with initial condition:
48 // T(x,y;1) = sin(pi*x) sin(pi*y)
49 // and boundary condition:
50 // T(boundary;t) = 0
51 
52 // which has the exact solution,
53 // z = exp(-K 2*pi^2 (1 - t)) * sin(pi*x) * sin(pi*y)
54 // which is the mirror image in time of the forward solution
55 
56 // For an adjoint consistent space-time formulation, the discrete
57 // adjoint can be obtained by marching backwards from the adjoint
58 // initial condition and solving the transpose of the discrete primal
59 // problem at the last nonlinear solve of the corresponding primal
60 // timestep. This necessitates the storage of the primal solution at
61 // all timesteps, which is accomplished here using either a
62 // MemorySolutionHistory or a FileSolutionHistory object. As the name suggests, this object
63 // simply stores the primal solution (and other vectors we may choose
64 // to save) in memory or disk, so that we can retrieve them later, whenever necessary.
65 
66 // The discrete adjoint system for implicit time steppers requires the
67 // localization of vectors other than system.solution, which is
68 // accomplished using the localize_vectors method. In this particular
69 // example, we use the localized adjoint solution to assemble the
70 // residual contribution for the current adjoint timestep from the last
71 // computed adjoint timestep.
72 
73 // Finally, The adjoint_advance_timestep method, the backwards time
74 // analog of advance_timestep prepares the time solver for solving the
75 // adjoint system, while the retrieve_timestep method retrieves the
76 // saved solutions at the current system.time, so that the adjoint
77 // sensitivity contribution for the current time can be computed.
78 
79 // Local includes
80 #include "initial.h"
81 #include "adjoint_initial.h"
82 #include "femparameters.h"
83 #include "heatsystem.h"
84 
85 // Libmesh includes
86 #include "libmesh/equation_systems.h"
87 #include "libmesh/dirichlet_boundaries.h"
88 #include "libmesh/system_norm.h"
89 #include "libmesh/numeric_vector.h"
90 #include "libmesh/enum_solver_package.h"
91 
92 #include "libmesh/mesh.h"
93 #include "libmesh/mesh_generation.h"
94 #include "libmesh/mesh_refinement.h"
95 
96 #include "libmesh/petsc_diff_solver.h"
97 #include "libmesh/steady_solver.h"
98 #include "libmesh/euler_solver.h"
99 #include "libmesh/euler2_solver.h"
100 #include "libmesh/newton_solver.h"
101 #include "libmesh/petsc_diff_solver.h"
102 #include "libmesh/twostep_time_solver.h"
103 
104 #include "libmesh/getpot.h"
105 #include "libmesh/tecplot_io.h"
106 #include "libmesh/gmv_io.h"
107 #include "libmesh/exodusII_io.h"
108 
109 #include "libmesh/sensitivity_data.h"
110 
111 // SolutionHistory Includes
112 #include "libmesh/solution_history.h"
113 #include "libmesh/memory_solution_history.h"
114 #include "libmesh/file_solution_history.h"
115 
116 // C++ includes
117 #include <iostream>
118 #include <sys/time.h>
119 #include <iomanip>
120 #include <memory>
121 
122 void write_output(EquationSystems & es,
123  unsigned int t_step, // The current time step count
124  std::string solution_type, // primal or adjoint solve
125  FEMParameters & param)
126 {
127  // Ignore parameters when there are no output formats available.
128  libmesh_ignore(es, t_step, solution_type, param);
129 
130 #ifdef LIBMESH_HAVE_GMV
131  if (param.output_gmv)
132  {
133  MeshBase & mesh = es.get_mesh();
134 
135  std::ostringstream file_name_gmv;
136  file_name_gmv << solution_type
137  << ".out.gmv."
138  << std::setw(2)
139  << std::setfill('0')
140  << std::right
141  << t_step;
142 
143  GMVIO(mesh).write_equation_systems(file_name_gmv.str(), es);
144  }
145 #endif
146 
147 #ifdef LIBMESH_HAVE_EXODUS_API
148  if (param.output_exodus)
149  {
150  MeshBase & mesh = es.get_mesh();
151 
152  // We write out one file per timestep. The files are named in
153  // the following way:
154  // foo.e
155  // foo.e-s002
156  // foo.e-s003
157  // ...
158  // so that, if you open the first one with Paraview, it actually
159  // opens the entire sequence of adapted files.
160  std::ostringstream file_name_exodus;
161 
162  file_name_exodus << solution_type << ".e";
163  if (t_step > 0)
164  file_name_exodus << "-s"
165  << std::setw(3)
166  << std::setfill('0')
167  << std::right
168  << t_step + 1;
169 
170  // TODO: Get the current time from the System...
171  ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
172  es,
173  1,
174  /*time=*/t_step + 1);
175  }
176 #endif
177 }
178 
180 {
181  // Use the prescribed FE type
182  system.fe_family() = param.fe_family[0];
183  system.fe_order() = param.fe_order[0];
184 
185  // Use analytical jacobians?
186  system.analytic_jacobians() = param.analytic_jacobians;
187 
188  // Verify analytic jacobians against numerical ones?
191 
192  // More desperate debugging options
194  system.print_solutions = param.print_solutions;
196  system.print_residuals = param.print_residuals;
198  system.print_jacobians = param.print_jacobians;
201 
202  // Solve this as a time-dependent or steady system
203  if (param.transient)
204  {
205  std::unique_ptr<UnsteadySolver> innersolver;
206  if (param.timesolver_core == "euler2")
207  {
208  auto euler2solver = std::make_unique<Euler2Solver>(system);
209  euler2solver->theta = param.timesolver_theta;
210  innersolver = std::move(euler2solver);
211  }
212  else if (param.timesolver_core == "euler")
213  {
214  auto eulersolver = std::make_unique<EulerSolver>(system);
215  eulersolver->theta = param.timesolver_theta;
216  innersolver = std::move(eulersolver);
217  }
218  else
219  libmesh_error_msg("This example (and unsteady adjoints in libMesh) only support Backward Euler and explicit methods.");
220 
221  if (param.timesolver_tolerance)
222  {
223  system.time_solver = std::make_unique<TwostepTimeSolver>(system);
224  auto timesolver = cast_ptr<TwostepTimeSolver *>(system.time_solver.get());
225 
226  timesolver->max_growth = param.timesolver_maxgrowth;
227  timesolver->target_tolerance = param.timesolver_tolerance;
228  timesolver->upper_tolerance = param.timesolver_upper_tolerance;
229  timesolver->component_norm = SystemNorm(param.timesolver_norm);
230  timesolver->core_time_solver = std::move(innersolver);
231 
232  // The Memory/File Solution History object we will set the system SolutionHistory object to
233  if(param.solution_history_type == "memory")
234  {
235  MemorySolutionHistory heatsystem_solution_history(system);
236  system.time_solver->set_solution_history(heatsystem_solution_history);
237  }
238  else if (param.solution_history_type == "file")
239  {
240  FileSolutionHistory heatsystem_solution_history(system);
241  system.time_solver->set_solution_history(heatsystem_solution_history);
242  }
243  else
244  libmesh_error_msg("Unrecognized solution history type: " << param.solution_history_type);
245  }
246  else
247  {
248  system.time_solver = std::move(innersolver);
249 
250  // The Memory/File Solution History object we will set the system SolutionHistory object to
251  if(param.solution_history_type == "memory")
252  {
253  MemorySolutionHistory heatsystem_solution_history(system);
254  system.time_solver->set_solution_history(heatsystem_solution_history);
255  }
256  else if (param.solution_history_type == "file")
257  {
258  FileSolutionHistory heatsystem_solution_history(system);
259  system.time_solver->set_solution_history(heatsystem_solution_history);
260  }
261  else
262  libmesh_error_msg("Unrecognized solution history type: " << param.solution_history_type);
263  }
264  }
265  else // !param.transient
266  system.time_solver = std::make_unique<SteadySolver>(system);
267 
268  system.time_solver->reduce_deltat_on_diffsolver_failure = param.deltat_reductions;
269  system.time_solver->quiet = param.time_solver_quiet;
270 
271  #ifdef LIBMESH_ENABLE_DIRICHLET
272  // Create any Dirichlet boundary conditions
273  typedef
274  std::map<boundary_id_type, FunctionBase<Number> *>::
275  const_iterator Iter;
276 
277  for (Iter i = param.dirichlet_conditions.begin();
278  i != param.dirichlet_conditions.end(); ++i)
279  {
280  boundary_id_type b = i->first;
281  FunctionBase<Number> *f = i->second;
282 
283  system.get_dof_map().add_dirichlet_boundary(DirichletBoundary({b},
285  f));
286 
287  libMesh::out << "Added Dirichlet boundary " << b << " for variables ";
288  for (std::size_t vi=0; vi != param.dirichlet_condition_variables[b].size(); ++vi)
290  libMesh::out << std::endl;
291  }
292  #endif // LIBMESH_ENABLE_DIRICHLET
293 
294  // Set the time stepping options
295  system.deltat = param.deltat;
296 
297  // And the integration options
299 
300  // And the nonlinear solver options
301  if (param.use_petsc_snes)
302  {
303 #ifdef LIBMESH_HAVE_PETSC
304  system.time_solver->diff_solver() = std::make_unique<PetscDiffSolver>(system);
305 #else
306  libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
307 #endif
308  }
309  else
310  {
311  system.time_solver->diff_solver() = std::make_unique<NewtonSolver>(system);
312  auto solver = cast_ptr<NewtonSolver *>(system.time_solver->diff_solver().get());
313 
314  solver->quiet = param.solver_quiet;
315  solver->verbose = param.solver_verbose;
316  solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
317  solver->minsteplength = param.min_step_length;
318  solver->relative_step_tolerance = param.relative_step_tolerance;
319  solver->relative_residual_tolerance = param.relative_residual_tolerance;
320  solver->require_residual_reduction = param.require_residual_reduction;
321  solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
322  if (system.time_solver->reduce_deltat_on_diffsolver_failure)
323  {
324  solver->continue_after_max_iterations = true;
325  solver->continue_after_backtrack_failure = true;
326  }
328 
329  // And the linear solver options
330  solver->max_linear_iterations = param.max_linear_iterations;
331  solver->initial_linear_tolerance = param.initial_linear_tolerance;
332  solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
333  }
334 }
335 
336 // The main program.
337 int main (int argc, char ** argv)
338 {
339  // Skip adaptive examples on a non-adaptive libMesh build
340 #ifndef LIBMESH_ENABLE_AMR
341  libmesh_ignore(argc, argv);
342  libmesh_example_requires(false, "--enable-amr");
343 #else
344  // Skip this 2D example if libMesh was compiled as 1D-only.
345  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
346 
347  // We use Dirichlet boundary conditions here
348 #ifndef LIBMESH_ENABLE_DIRICHLET
349  libmesh_example_requires(false, "--enable-dirichlet");
350 #endif
351 
352  // Initialize libMesh.
353  LibMeshInit init (argc, argv);
354 
355  // This doesn't converge with Trilinos for some reason...
356  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
357 
358  libMesh::out << "Started " << argv[0] << std::endl;
359 
360  // Make sure the general input file exists, and parse it
361  {
362  std::ifstream i("general.in");
363  libmesh_error_msg_if(!i, '[' << init.comm().rank() << "] Can't find general.in; exiting early.");
364  }
365  GetPot infile("general.in");
366 
367  // But allow the command line to override it.
368  infile.parse_command_line(argc, argv);
369 
370  // Read in parameters from the input file
371  FEMParameters param(init.comm());
372  param.read(infile);
373 
374  // Create a mesh with the given dimension, distributed
375  // across the default MPI communicator.
376  Mesh mesh(init.comm(), cast_int<unsigned char>(param.dimension));
377 
378  // And an object to refine it
379  auto mesh_refinement = std::make_unique<MeshRefinement>(mesh);
380 
381  // And an EquationSystems to run on it
382  EquationSystems equation_systems (mesh);
383 
384  libMesh::out << "Building mesh" << std::endl;
385 
386  // Build a unit square
387  ElemType elemtype;
388 
389  if (param.elementtype == "tri" ||
390  param.elementtype == "unstructured")
391  elemtype = TRI3;
392  else
393  elemtype = QUAD4;
394 
395  MeshTools::Generation::build_square (mesh, param.coarsegridx, param.coarsegridy,
396  param.domain_xmin, param.domain_xmin + param.domain_edge_width,
397  param.domain_ymin, param.domain_ymin + param.domain_edge_length,
398  elemtype);
399 
400  libMesh::out << "Building system" << std::endl;
401 
402  HeatSystem & system = equation_systems.add_system<HeatSystem> ("HeatSystem");
403 
404  set_system_parameters(system, param);
405 
406  libMesh::out << "Initializing systems" << std::endl;
407 
408  // Initialize the system
409  equation_systems.init ();
410 
411  // Refine the grid again if requested
412  for (unsigned int i=0; i != param.extrarefinements; ++i)
413  {
414  mesh_refinement->uniformly_refine(1);
415  equation_systems.reinit();
416  }
417 
418  libMesh::out << "Setting primal initial conditions" << std::endl;
419 
421 
423  equation_systems.parameters);
424 
425  // Output the H1 norm of the initial conditions
426  libMesh::out << "|U("
427  << system.time
428  << ")|= "
429  << system.calculate_norm(*system.solution, 0, H1)
430  << std::endl
431  << std::endl;
432 
433  // Tell system that we have two qois
434  system.init_qois(1);
435 
436  // Add the adjoint vectors (and the old adjoint vectors) to system
437  system.time_solver->init_adjoints();
438 
439  // Close up any resources initial.C needed
441 
442  // Plot the initial conditions
443  write_output(equation_systems, 0, "primal", param);
444 
445  // Print information about the mesh and system to the screen.
446  mesh.print_info();
447  equation_systems.print_info();
448 
449  // In optimized mode we catch any solver errors, so that we can
450  // write the proper footers before closing. In debug mode we just
451  // let the exception throw so that gdb can grab it.
452 #if defined(NDEBUG) && defined(LIBMESH_ENABLE_EXCEPTIONS)
453  try
454 #endif
455  {
456  // Now we begin the timestep loop to compute the time-accurate
457  // solution of the equations.
458  for (unsigned int t_step=param.initial_timestep;
459  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
460  {
461  // A pretty update message
462  libMesh::out << " Solving time step "
463  << t_step
464  << ", time = "
465  << system.time
466  << std::endl;
467 
468  // Solve the forward problem at time t, to obtain the solution at time t + dt
469  system.solve();
470 
471 
472  // Output the H1 norm of the computed solution
473  libMesh::out << "|U("
474  << system.time + system.time_solver->last_completed_timestep_size()
475  << ")|= "
476  << system.calculate_norm(*system.solution, 0, H1)
477  << std::endl;
478 
479  // Advance to the next timestep in a transient problem
480  libMesh::out << "Advancing timestep" << std::endl << std::endl;
481  system.time_solver->advance_timestep();
482 
483  // Write out this timestep
484  write_output(equation_systems, t_step+1, "primal", param);
485  }
486  // End timestep loop
487 
489 
490  // Now we will solve the backwards in time adjoint problem
491  libMesh::out << std::endl << "Solving the adjoint problem" << std::endl;
492 
493  // We need to tell the library that it needs to project the adjoint, so
494  // MemorySolutionHistory knows it has to save it
495 
496  const std::string & adjoint_solution_name0 = "adjoint_solution0";
497  const std::string & old_adjoint_solution_name0 = "_old_adjoint_solution0";
498  // Tell the library to project the adjoint vector, and hence, memory solution history to
499  // save it
500  system.set_vector_preservation(adjoint_solution_name0, true);
501  system.set_vector_preservation(old_adjoint_solution_name0, true);
502 
503  libMesh::out << "Setting adjoint initial conditions Z("
504  << system.time
505  << ")"
506  <<std::endl;
507 
508  // The first thing we have to do is to apply the adjoint initial
509  // condition. The user should supply these. Here they are specified
510  // in the functions adjoint_initial_value and adjoint_initial_gradient
513  equation_systems.parameters,
514  system.get_adjoint_solution(0));
515 
516  // Since we have specified an adjoint solution for the current
517  // time (T), set the adjoint_already_solved boolean to true, so
518  // we dont solve unnecessarily in the adjoint sensitivity method
519  system.set_adjoint_already_solved(true);
520 
521  Real Z_norm = system.calculate_norm(system.get_adjoint_solution(), 0, H1);
522 
523  libMesh::out << "|Z("
524  << system.time
525  << ")|= "
526  << Z_norm
527  << std::endl
528  << std::endl;
529 
530  // Call the adjoint advance timestep here to save the adjoint
531  // initial conditions to disk
532  system.time_solver->adjoint_advance_timestep();
533 
534  Real Z_old_norm = system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1);
535 
536  // Make sure that we have copied the old adjoint solution properly
537  libmesh_assert(Z_norm == Z_old_norm);
538 
539  libMesh::out << "|Z_old("
540  << system.time + (system.time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
541  << ")|= "
542  << Z_old_norm
543  << std::endl
544  << std::endl;
545 
546  write_output(equation_systems, param.n_timesteps, "dual", param);
547 
548  // Now that the adjoint initial condition is set, we will start the
549  // backwards in time adjoint integration
550 
551  // For loop stepping backwards in time
552  for (unsigned int t_step=param.initial_timestep;
553  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
554  {
555  //A pretty update message
556  libMesh::out << " Solving adjoint time step "
557  << t_step
558  << ", time = "
559  << system.time
560  << std::endl;
561 
562  // Output the H1 norm of the retrieved primal solution from the last call
563  // to adjoint_advance_timestep
564  libMesh::out << "|U("
565  << system.time
566  << ")|= "
567  << system.calculate_norm(*system.solution, 0, H1)
568  << std::endl;
569 
570  libMesh::out << "|U("
571  << system.time - (system.time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
572  << ")|= "
573  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
574  << std::endl;
575 
576  system.set_adjoint_already_solved(false);
577 
578  system.adjoint_solve();
579 
580  Z_norm = system.calculate_norm(system.get_adjoint_solution(), 0, H1);
581 
582  libMesh::out << "|Z("
583  << system.time
584  << ")|= "
585  << Z_norm
586  << std::endl
587  << std::endl;
588 
589  // Now that we have solved the adjoint, set the
590  // adjoint_already_solved boolean to true, so we dont solve
591  // unnecessarily in the error estimator
592  system.set_adjoint_already_solved(true);
593 
594  libMesh::out << "Saving adjoint and retrieving primal solutions at time t=" << system.time - system.deltat << std::endl;
595 
596  // The adjoint_advance_timestep function calls the retrieve and store
597  // function of the memory_solution_history class via the
598  // memory_solution_history object we declared earlier. The
599  // retrieve function sets the system primal vectors to their
600  // values at the current time instance.
601  system.time_solver->adjoint_advance_timestep();
602 
603  Z_old_norm = system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1);
604 
605  // Make sure that we have copied the old adjoint solution properly
606  libmesh_assert(Z_norm == Z_old_norm);
607 
608  libMesh::out << "|Z_old("
609  << system.time + (system.time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
610  << ")|= "
611  << Z_old_norm
612  << std::endl
613  << std::endl;
614 
615  // Get a pointer to the primal solution vector
616  NumericVector<Number> & primal_solution = *system.solution;
617 
618  // Get a pointer to the solution vector of the adjoint problem for QoI 0
619  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
620 
621  // Swap the primal and dual solutions so we can write out the adjoint solution
622  primal_solution.swap(dual_solution_0);
623 
624  write_output(equation_systems, param.n_timesteps - (t_step + 1), "dual", param);
625 
626  // Swap back
627  primal_solution.swap(dual_solution_0);
628  }
629  // End adjoint timestep loop
630 
631  // Now that we have computed both the primal and adjoint solutions, we compute the sensitivities to the parameter p
632  // dQ/dp = int_{0}^{T} partialQ/partialp - partialR/partialp(u,z;p) dt
633  // The quantity partialQ/partialp - partialR/partialp(u,z;p) is evaluated internally by the ImplicitSystem::adjoint_qoi_parameter_sensitivity function.
634  // This sensitivity evaluation is called internally by an overloaded TimeSolver::integrate_adjoint_sensitivity method which we call below.
635 
636  // Reset the time
637  system.time = 0.0;
638 
639  // Prepare the quantities we need to pass to TimeSolver::integrate_adjoint_sensitivity
640  QoISet qois;
641 
642  qois.add_indices({0});
643  qois.set_weight(0, 1.0);
644 
645  // A SensitivityData object
646  SensitivityData sensitivities(qois, system, system.get_parameter_vector());
647 
648  // Accumulator for the time integrated total sensitivity
649  Number total_sensitivity = 0.0;
650 
651  // Retrieve the primal and adjoint solutions at the current timestep
652  system.time_solver->retrieve_timestep();
653 
654  Z_old_norm = system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1);
655 
656  // Assert that the old adjoint values match hard coded numbers for 1st timestep or 1st half time step from the adjoint solve
657  if(param.timesolver_tolerance)
658  {
659  libmesh_error_msg_if(std::abs(Z_old_norm - (2.23366)) >= 2.e-4,
660  "Mismatch in expected Z0_old norm for the 1st half timestep!");
661  }
662  else
663  {
664  libmesh_error_msg_if(std::abs(Z_old_norm - (2.23627)) >= 2.e-4,
665  "Mismatch in expected Z0_old norm for the 1st timestep!");
666  }
667 
668  // A pretty update message
669  libMesh::out << "Retrieved, "
670  << "time = "
671  << system.time
672  << std::endl;
673 
674  libMesh::out << "|U("
675  << system.time
676  << ")|= "
677  << system.calculate_norm(*system.solution, 0, H1)
678  << std::endl;
679 
680  libMesh::out << "|U_old("
681  << system.time
682  << ")|= "
683  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
684  << std::endl;
685 
686  libMesh::out << "|Z("
687  << system.time
688  << ")|= "
689  << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
690  << std::endl
691  << std::endl;
692 
693  libMesh::out << "|Z_old("
694  << system.time + (system.time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
695  << ")|= "
696  << system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1)
697  << std::endl
698  << std::endl;
699 
700  // Now we begin the timestep loop to compute the time-accurate
701  // adjoint sensitivities
702  for (unsigned int t_step=param.initial_timestep;
703  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
704  {
705  // Call the postprocess function which we have overloaded to compute
706  // accumulate the perturbed residuals
707  system.time_solver->integrate_adjoint_sensitivity(qois, system.get_parameter_vector(), sensitivities);
708 
709  // A pretty update message
710  libMesh::out << "Retrieved, "
711  << "time = "
712  << system.time
713  << std::endl;
714 
715  libMesh::out << "|U("
716  << system.time
717  << ")|= "
718  << system.calculate_norm(*system.solution, 0, H1)
719  << std::endl;
720 
721  libMesh::out << "|U_old("
722  << system.time
723  << ")|= "
724  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
725  << std::endl;
726 
727  libMesh::out << "|Z("
728  << system.time
729  << ")|= "
730  << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
731  << std::endl
732  << std::endl;
733 
734  libMesh::out << "|Z_old("
735  << system.time
736  << ")|= "
737  << system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1)
738  << std::endl
739  << std::endl;
740 
741  // Add the contribution of the current timestep to the total sensitivity
742  total_sensitivity += sensitivities[0][0];
743  }
744 
745  // Print it out
746  libMesh::out << "Sensitivity of QoI 0 w.r.t parameter 0 is: "
747  << total_sensitivity
748  << std::endl;
749 
750  // Hard coded test to ensure that the actual numbers we are
751  // getting are what they should be
752  // The 2e-4 tolerance is chosen to ensure success even with
753  // 32-bit floats
754  if(param.timesolver_tolerance)
755  {
756  libmesh_error_msg_if(std::abs(system.time - (1.0089)) >= 2.e-4,
757  "Mismatch in end time reached by adaptive timestepper!");
758 
759  libmesh_error_msg_if(std::abs(total_sensitivity - 4.87767) >= 3.e-3,
760  "Mismatch in sensitivity gold value!");
761  }
762  else
763  {
764  libmesh_error_msg_if(std::abs(total_sensitivity - 4.83551) >= 2.e-4,
765  "Mismatch in sensitivity gold value!");
766  }
767  }
768 #if defined(NDEBUG) && defined(LIBMESH_ENABLE_EXCEPTIONS)
769  catch (...)
770  {
771  libMesh::err << '[' << mesh.processor_id()
772  << "] Caught exception; exiting early." << std::endl;
773  }
774 #endif
775 
776  libMesh::err << '[' << mesh.processor_id()
777  << "] Completing output."
778  << std::endl;
779 
780  // All done.
781  return 0;
782 
783 #endif // LIBMESH_ENABLE_AMR
784 }
OStreamProxy err
libMesh::Real timesolver_upper_tolerance
Definition: femparameters.h:39
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1595
bool print_solution_norms
Number adjoint_initial_value(const Point &p, const Parameters &, const std::string &, const std::string &)
ElemType
Defines an enum for geometric element types.
libMesh::Real deltat
Definition: femparameters.h:39
void write_output(EquationSystems &es, unsigned int t_step, std::string solution_type, FEMParameters &param)
Definition: adjoints_ex5.C:122
bool print_element_residuals
std::string solution_history_type
Definition: femparameters.h:38
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:215
libMesh::Real timesolver_tolerance
Definition: femparameters.h:39
void finish_initialization()
Definition: initial.C:11
void set_adjoint_already_solved(bool setting)
Setter for the adjoint_already_solved boolean.
Definition: system.h:412
bool analytic_jacobians
virtual void set_constrain_in_solver(bool enable)
set_constrain_in_solver to false to apply constraints only via residual terms in the systems to be so...
Definition: diff_system.C:422
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Definition: diff_system.h:359
Gradient adjoint_initial_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
int extra_quadrature_order
A member int that can be employed to indicate increased or reduced quadrature order.
Definition: system.h:1558
Real numerical_jacobian_h
If calculating numeric jacobians is required, the FEMSystem will perturb each solution vector entry b...
Definition: fem_system.h:187
std::string timesolver_core
Definition: femparameters.h:37
bool require_residual_reduction
std::unique_ptr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:251
MeshBase & mesh
bool print_residual_norms
void init_qois(unsigned int n_qois)
Accessors for qoi and qoi_error_estimates vectors.
Definition: system.C:2319
libMesh::Real timesolver_maxgrowth
Definition: femparameters.h:39
int main(int argc, char **argv)
Definition: adjoints_ex5.C:337
libMesh::Real relative_step_tolerance
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:364
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:374
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
SolverPackage default_solver_package()
Definition: libmesh.C:1050
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr) const
Projects arbitrary functions onto the current solution.
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:338
void libmesh_ignore(const Args &...)
unsigned int max_linear_iterations
int8_t boundary_id_type
Definition: id_types.h:51
bool constrain_in_solver
bool print_residual_norms
Set print_residual_norms to true to print |F| whenever it is assembled.
Definition: diff_system.h:349
double minimum_linear_tolerance
ParameterVector & get_parameter_vector()
Definition: heatsystem.h:76
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:278
std::vector< std::string > fe_family
libMesh::Real verify_analytic_jacobians
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1573
Real calculate_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type, std::set< unsigned int > *skip_dimensions=nullptr) const
Definition: system.C:1672
bool print_residuals
Set print_residuals to true to print F whenever it is assembled.
Definition: diff_system.h:354
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
libmesh_assert(ctx)
double initial_linear_tolerance
int extra_quadrature_order
void read(GetPot &input, const std::vector< std::string > *other_variable_names=nullptr)
bool print_solutions
Set print_solutions to true to print U whenever it is used in an assembly() call. ...
Definition: diff_system.h:344
libMesh::Real relative_residual_tolerance
void set_system_parameters(HeatSystem &system, FEMParameters &param)
Definition: adjoints_ex5.C:179
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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:1087
bool print_jacobian_norms
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
OStreamProxy out
libMesh::Real min_step_length
std::string & fe_family()
Definition: heatsystem.h:46
libMesh::Real timesolver_theta
Definition: femparameters.h:39
Gradient initial_grad(const Point &, const Parameters &, const std::string &, const std::string &)
Definition: initial.C:30
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:1193
bool print_element_jacobians
Set print_element_jacobians to true to print each J_elem contribution.
Definition: diff_system.h:379
bool print_element_jacobians
unsigned int max_nonlinear_iterations
double linear_tolerance_multiplier
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:94
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr, int is_adjoint=-1) const
Projects arbitrary functions onto a vector of degree of freedom values for the current system...
virtual void solve() override
Invokes the solver associated with the system.
Definition: fem_system.C:1076
bool time_solver_quiet
bool & analytic_jacobians()
Definition: heatsystem.h:49
void read_initial_parameters()
Definition: initial.C:7
libMesh::Real numerical_jacobian_h
std::map< libMesh::boundary_id_type, libMesh::FunctionBase< libMesh::Number > * > dirichlet_conditions
Definition: femparameters.h:91
const DofMap & get_dof_map() const
Definition: system.h:2293
std::vector< libMesh::FEMNormType > timesolver_norm
Definition: femparameters.h:42
template class LIBMESH_EXPORT NumericVector< Number >
Number initial_value(const Point &, const Parameters &, const std::string &, const std::string &)
Definition: initial.C:18
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) override
This function sets the _is_adjoint boolean member of TimeSolver to true and then calls the adjoint_so...
Definition: diff_system.C:150
unsigned int & fe_order()
Definition: heatsystem.h:47
const NumericVector< Number > & get_vector(std::string_view vec_name) const
Definition: system.C:918
unsigned int deltat_reductions
Definition: femparameters.h:36