libMesh
adjoints_ex2.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 // <h1>Adjoints Example 2 - Laplace Equation in the L-Shaped Domain with
21 // Adjoint based sensitivity</h1>
22 // \author Roy Stogner
23 // \date 2003
24 //
25 // This example solves the Laplace equation on the classic "L-shaped"
26 // domain with adaptive mesh refinement. The exact solution is
27 // u(r,\theta) = r^{2/3} * \sin ( (2/3) * \theta). We scale this
28 // exact solution by a combination of parameters, (\alpha_{1} + 2
29 // *\alpha_{2}) to get u = (\alpha_{1} + 2 *\alpha_{2}) * r^{2/3} *
30 // \sin ( (2/3) * \theta), which again satisfies the Laplace
31 // Equation. We define the Quantity of Interest in element_qoi.C, and
32 // compute the sensitivity of the QoI to \alpha_{1} and \alpha_{2}
33 // using the adjoint sensitivity method. Since we use the adjoint
34 // capabilities of libMesh in this example, we use the DiffSystem
35 // framework. This file (main.C) contains the declaration of mesh and
36 // equation system objects, L-shaped.C contains the assembly of the
37 // system, element_qoi_derivative.C contains
38 // the RHS for the adjoint system. Postprocessing to compute the the
39 // QoIs is done in element_qoi.C
40 
41 // The initial mesh contains three QUAD9 elements which represent the
42 // standard quadrants I, II, and III of the domain [-1,1]x[-1,1],
43 // i.e.
44 // Element 0: [-1,0]x[ 0,1]
45 // Element 1: [ 0,1]x[ 0,1]
46 // Element 2: [-1,0]x[-1,0]
47 // The mesh is provided in the standard libMesh ASCII format file
48 // named "lshaped.xda". In addition, an input file named "general.in"
49 // is provided which allows the user to set several parameters for
50 // the solution so that the problem can be re-run without a
51 // re-compile. The solution technique employed is to have a
52 // refinement loop with a linear (forward and adjoint) solve inside followed by a
53 // refinement of the grid and projection of the solution to the new grid
54 // In the final loop iteration, there is no additional
55 // refinement after the solve. In the input file "general.in", the variable
56 // "max_adaptivesteps" controls the number of refinement steps, and
57 // "refine_fraction" / "coarsen_fraction" determine the number of
58 // elements which will be refined / coarsened at each step.
59 
60 // C++ includes
61 #include <iostream>
62 #include <iomanip>
63 #include <memory>
64 
65 // General libMesh includes
66 #include "libmesh/eigen_sparse_linear_solver.h"
67 #include "libmesh/enum_solver_package.h"
68 #include "libmesh/enum_solver_type.h"
69 #include "libmesh/equation_systems.h"
70 #include "libmesh/error_vector.h"
71 #include "libmesh/mesh.h"
72 #include "libmesh/mesh_refinement.h"
73 #include "libmesh/newton_solver.h"
74 #include "libmesh/numeric_vector.h"
75 #include "libmesh/petsc_diff_solver.h"
76 #include "libmesh/steady_solver.h"
77 #include "libmesh/system_norm.h"
78 
79 // Sensitivity Calculation related includes
80 #include "libmesh/parameter_vector.h"
81 #include "libmesh/sensitivity_data.h"
82 
83 // Error Estimator includes
84 #include "libmesh/kelly_error_estimator.h"
85 #include "libmesh/patch_recovery_error_estimator.h"
86 
87 // Adjoint Related includes
88 #include "libmesh/adjoint_residual_error_estimator.h"
89 #include "libmesh/qoi_set.h"
90 
91 // libMesh I/O includes
92 #include "libmesh/getpot.h"
93 #include "libmesh/gmv_io.h"
94 #include "libmesh/exodusII_io.h"
95 
96 // Local includes
97 #include "femparameters.h"
98 #include "L-shaped.h"
99 #include "L-qoi.h"
100 
101 // Bring in everything from the libMesh namespace
102 using namespace libMesh;
103 
104 // Local function declarations
105 
106 // Number output files, the files are give a prefix of primal or adjoint_i depending on
107 // whether the output is the primal solution or the dual solution for the ith QoI
108 
109 // Write gmv output
111  unsigned int a_step, // The adaptive step count
112  std::string solution_type, // primal or adjoint solve
113  FEMParameters & param)
114 {
115  // Ignore parameters when there are no output formats available.
116  libmesh_ignore(es, a_step, solution_type, param);
117 
118 #ifdef LIBMESH_HAVE_GMV
119  if (param.output_gmv)
120  {
121  MeshBase & mesh = es.get_mesh();
122 
123  std::ostringstream file_name_gmv;
124  file_name_gmv << solution_type
125  << ".out.gmv."
126  << std::setw(2)
127  << std::setfill('0')
128  << std::right
129  << a_step;
130 
131  GMVIO(mesh).write_equation_systems(file_name_gmv.str(), es);
132  }
133 #endif
134 
135 #ifdef LIBMESH_HAVE_EXODUS_API
136  if (param.output_exodus)
137  {
138  MeshBase & mesh = es.get_mesh();
139 
140  // We write out one file per adaptive step. The files are named in
141  // the following way:
142  // foo.e
143  // foo.e-s002
144  // foo.e-s003
145  // ...
146  // so that, if you open the first one with Paraview, it actually
147  // opens the entire sequence of adapted files.
148  std::ostringstream file_name_exodus;
149 
150  file_name_exodus << solution_type << ".e";
151  if (a_step > 0)
152  file_name_exodus << "-s"
153  << std::setw(3)
154  << std::setfill('0')
155  << std::right
156  << a_step + 1;
157 
158  // We write each adaptive step as a pseudo "time" step, where the
159  // time simply matches the (1-based) adaptive step we are on.
160  ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
161  es,
162  1,
163  /*time=*/a_step + 1);
164  }
165 #endif
166 }
167 
168 
170 {
171  // Eigen's BiCGSTAB doesn't seem reliable at the full refinement
172  // level we use here.
173 #ifdef LIBMESH_HAVE_EIGEN_SPARSE
174  EigenSparseLinearSolver<Number> * eigen_linear_solver =
175  dynamic_cast<EigenSparseLinearSolver<Number> *>(&linear_solver);
176 
177  if (eigen_linear_solver)
178  eigen_linear_solver->set_solver_type(SPARSELU);
179 #else
180  libmesh_ignore(linear_solver);
181 #endif
182 }
183 
185 {
186  auto diff_solver = cast_ptr<NewtonSolver*>(system.get_time_solver().diff_solver().get());
187  if (diff_solver) // Some compilers don't like dynamic cast of nullptr?
188  {
189  auto solver = cast_ptr<NewtonSolver*>(diff_solver);
190  if (solver)
191  adjust_linear_solver(solver->get_linear_solver());
192  }
193 
194  LinearSolver<Number> * linear_solver = system.get_linear_solver();
195  if (linear_solver)
196  adjust_linear_solver(*linear_solver);
197 }
198 
199 
200 // Set the parameters for the nonlinear and linear solvers to be used during the simulation
202 {
203  // Use analytical jacobians?
204  system.analytic_jacobians() = param.analytic_jacobians;
205 
206  // Verify analytic jacobians against numerical ones?
208 
209  // Use the prescribed FE type
210  system.fe_family() = param.fe_family[0];
211  system.fe_order() = param.fe_order[0];
212 
213  // More desperate debugging options
215  system.print_solutions = param.print_solutions;
217  system.print_residuals = param.print_residuals;
219  system.print_jacobians = param.print_jacobians;
220 
221  // No transient time solver
222  system.time_solver = std::make_unique<SteadySolver>(system);
223 
224  // Nonlinear solver options
225  if (param.use_petsc_snes)
226  {
227 #ifdef LIBMESH_HAVE_PETSC
228  system.time_solver->diff_solver() = std::make_unique<PetscDiffSolver>(system);
229 #else
230  libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
231 #endif
232  }
233  else
234  {
235  system.time_solver->diff_solver() = std::make_unique<NewtonSolver>(system);
236  auto solver = cast_ptr<NewtonSolver*>(system.time_solver->diff_solver().get());
237 
238  solver->quiet = param.solver_quiet;
239  solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
240  solver->minsteplength = param.min_step_length;
241  solver->relative_step_tolerance = param.relative_step_tolerance;
242  solver->relative_residual_tolerance = param.relative_residual_tolerance;
243  solver->require_residual_reduction = param.require_residual_reduction;
244  solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
245  if (system.time_solver->reduce_deltat_on_diffsolver_failure)
246  {
247  solver->continue_after_max_iterations = true;
248  solver->continue_after_backtrack_failure = true;
249  }
251 
252  // And the linear solver options
253  solver->max_linear_iterations = param.max_linear_iterations;
254  solver->initial_linear_tolerance = param.initial_linear_tolerance;
255  solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
256  adjust_linear_solvers(system);
257  }
258 }
259 
260 // Build the mesh refinement object and set parameters for refining/coarsening etc
261 
262 #ifdef LIBMESH_ENABLE_AMR
263 
264 std::unique_ptr<MeshRefinement> build_mesh_refinement(MeshBase & mesh,
265  FEMParameters & param)
266 {
267  auto mesh_refinement = std::make_unique<MeshRefinement>(mesh);
268  mesh_refinement->coarsen_by_parents() = true;
269  mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
270  mesh_refinement->nelem_target() = param.nelem_target;
271  mesh_refinement->refine_fraction() = param.refine_fraction;
272  mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
273  mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
274 
275  return mesh_refinement;
276 }
277 
278 #endif // LIBMESH_ENABLE_AMR
279 
280 // This is where we declare the error estimators to be built and used for
281 // mesh refinement. The adjoint residual estimator needs two estimators.
282 // One for the forward component of the estimate and one for the adjoint
283 // weighting factor. Here we use the Patch Recovery indicator to estimate both the
284 // forward and adjoint weights. The H1 seminorm component of the error is used
285 // as dictated by the weak form the Laplace equation.
286 
287 std::unique_ptr<ErrorEstimator> build_error_estimator(FEMParameters & param)
288 {
289  if (param.indicator_type == "kelly")
290  {
291  libMesh::out << "Using Kelly Error Estimator" << std::endl;
292 
293  return std::make_unique<KellyErrorEstimator>();
294  }
295  else if (param.indicator_type == "adjoint_residual")
296  {
297  libMesh::out << "Using Adjoint Residual Error Estimator with Patch Recovery Weights" << std::endl << std::endl;
298 
299  auto adjoint_residual_estimator = std::make_unique<AdjointResidualErrorEstimator>();
300 
301  adjoint_residual_estimator->error_plot_suffix = "error.gmv";
302 
303  adjoint_residual_estimator->primal_error_estimator() = std::make_unique<PatchRecoveryErrorEstimator>();
304  adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
305 
306  adjoint_residual_estimator->dual_error_estimator() = std::make_unique<PatchRecoveryErrorEstimator>();
307  adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
308 
309  return adjoint_residual_estimator;
310  }
311  else
312  libmesh_error_msg("Unknown indicator_type = " << param.indicator_type);
313 }
314 
315 // The main program.
316 int main (int argc, char ** argv)
317 {
318  // Initialize libMesh.
319  LibMeshInit init (argc, argv);
320 
321  // This example requires a linear solver package.
322  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
323  "--enable-petsc, --enable-trilinos, or --enable-eigen");
324 
325  // Skip adaptive examples on a non-adaptive libMesh build
326 #ifndef LIBMESH_ENABLE_AMR
327  libmesh_example_requires(false, "--enable-amr");
328 #else
329 
330  libMesh::out << "Started " << argv[0] << std::endl;
331 
332  // Make sure the general input file exists, and parse it
333  {
334  std::ifstream i("general.in");
335  libmesh_error_msg_if(!i, '[' << init.comm().rank() << "] Can't find general.in; exiting early.");
336  }
337 
338  // Read in parameters from the input file
339  GetPot infile("general.in");
340  // But allow the command line to override it.
341  infile.parse_command_line(argc, argv);
342 
343  FEMParameters param(init.comm());
344  param.read(infile);
345 
346  // Skip this default-2D example if libMesh was compiled as 1D-only.
347  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
348 
349  // Create a mesh, with dimension to be overridden later, distributed
350  // across the default MPI communicator.
351  Mesh mesh(init.comm());
352 
353  // And an object to refine it
354  std::unique_ptr<MeshRefinement> mesh_refinement =
355  build_mesh_refinement(mesh, param);
356 
357  // And an EquationSystems to run on it
358  EquationSystems equation_systems (mesh);
359 
360  libMesh::out << "Reading in and building the mesh" << std::endl;
361 
362  // Read in the mesh
363  mesh.read(param.domainfile.c_str());
364  // Make all the elements of the mesh second order so we can compute
365  // with a higher order basis
367 
368  // Create a mesh refinement object to do the initial uniform refinements
369  // on the coarse grid read in from lshaped.xda
370  MeshRefinement initial_uniform_refinements(mesh);
371  initial_uniform_refinements.uniformly_refine(param.coarserefinements);
372 
373  libMesh::out << "Building system" << std::endl;
374 
375  // Build the FEMSystem
376  LaplaceSystem & system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
377 
378  QoISet qois;
379 
380  qois.add_indices({0});
381 
382  qois.set_weight(0, 0.5);
383 
384  // Put some scope here to test that the cloning is working right
385  {
386  LaplaceQoI qoi;
387  system.attach_qoi(&qoi);
388  }
389 
390  // Set its parameters
391  set_system_parameters(system, param);
392 
393  libMesh::out << "Initializing systems" << std::endl;
394 
395  equation_systems.init ();
396 
397  // Print information about the mesh and system to the screen.
398  mesh.print_info();
399  equation_systems.print_info();
400 
401  {
402  // Adaptively solve the timestep
403  unsigned int a_step = 0;
404  for (; a_step != param.max_adaptivesteps; ++a_step)
405  {
406  // We can't adapt to both a tolerance and a
407  // target mesh size
408  if (param.global_tolerance != 0.)
409  libmesh_assert_equal_to (param.nelem_target, 0);
410  // If we aren't adapting to a tolerance we need a
411  // target mesh size
412  else
413  libmesh_assert_greater (param.nelem_target, 0);
414 
415  // Solve the forward problem
416  system.solve();
417 
418  // Write out the computed primal solution
419  write_output(equation_systems, a_step, "primal", param);
420 
421  // Get a pointer to the primal solution vector
422  NumericVector<Number> & primal_solution = *system.solution;
423 
424  // A SensitivityData object to hold the qois and parameters
425  SensitivityData sensitivities(qois, system, system.get_parameter_vector());
426 
427  // Make sure we get the contributions to the adjoint RHS from the sides
428  system.assemble_qoi_sides = true;
429 
430  // Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity
431  // function, so we have to set the adjoint_already_solved boolean to false
432  system.set_adjoint_already_solved(false);
433 
434  // Compute the sensitivities
435  system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
436 
437  // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
438  system.set_adjoint_already_solved(true);
439 
440  GetPot infile_l_shaped("l-shaped.in");
441 
442  Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
443  Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
444  Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
445  Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
446 
447  libMesh::out << "Adaptive step "
448  << a_step
449  << ", we have "
450  << mesh.n_active_elem()
451  << " active elements and "
452  << equation_systems.n_active_dofs()
453  << " active dofs."
454  << std::endl;
455 
456  libMesh::out << "Sensitivity of QoI one to Parameter one is "
457  << sensitivity_QoI_0_0_computed
458  << std::endl;
459  libMesh::out << "Sensitivity of QoI one to Parameter two is "
460  << sensitivity_QoI_0_1_computed
461  << std::endl;
462 
463  libMesh::out << "The relative error in sensitivity QoI_0_0 is "
464  << std::setprecision(17)
465  << std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact) / std::abs(sensitivity_QoI_0_0_exact)
466  << std::endl;
467 
468  libMesh::out << "The relative error in sensitivity QoI_0_1 is "
469  << std::setprecision(17)
470  << std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact) / std::abs(sensitivity_QoI_0_1_exact)
471  << std::endl
472  << std::endl;
473 
474  // Get a pointer to the solution vector of the adjoint problem for QoI 0
475  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
476 
477  // Swap the primal and dual solutions so we can write out the adjoint solution
478  primal_solution.swap(dual_solution_0);
479  write_output(equation_systems, a_step, "adjoint_0", param);
480 
481  // Swap back
482  primal_solution.swap(dual_solution_0);
483 
484  // We have to refine either based on reaching an error tolerance or
485  // a number of elements target, which should be verified above
486  // Otherwise we flag elements by error tolerance or nelem target
487 
488  // Uniform refinement
489  if (param.refine_uniformly)
490  {
491  libMesh::out << "Refining Uniformly" << std::endl << std::endl;
492 
493  mesh_refinement->uniformly_refine(1);
494  }
495  // Adaptively refine based on reaching an error tolerance
496  else if (param.global_tolerance >= 0. && param.nelem_target == 0.)
497  {
498  // Now we construct the data structures for the mesh refinement process
499  ErrorVector error;
500 
501  // Build an error estimator object
502  std::unique_ptr<ErrorEstimator> error_estimator =
503  build_error_estimator(param);
504 
505  // Estimate the error in each element using the Adjoint Residual or Kelly error estimator
506  error_estimator->estimate_error(system, error);
507 
508  mesh_refinement->flag_elements_by_error_tolerance (error);
509 
510  mesh_refinement->refine_and_coarsen_elements();
511  }
512  // Adaptively refine based on reaching a target number of elements
513  else
514  {
515  // Now we construct the data structures for the mesh refinement process
516  ErrorVector error;
517 
518  // Build an error estimator object
519  std::unique_ptr<ErrorEstimator> error_estimator =
520  build_error_estimator(param);
521 
522  // Estimate the error in each element using the Adjoint Residual or Kelly error estimator
523  error_estimator->estimate_error(system, error);
524 
525  if (mesh.n_active_elem() >= param.nelem_target)
526  {
527  libMesh::out<<"We reached the target number of elements."<<std::endl <<std::endl;
528  break;
529  }
530 
531  mesh_refinement->flag_elements_by_nelem_target (error);
532 
533  mesh_refinement->refine_and_coarsen_elements();
534  }
535 
536  // Dont forget to reinit the system after each adaptive refinement !
537  equation_systems.reinit();
538 
539  // Fix up the linear solver options if that reinit just cleared it
540  adjust_linear_solvers(system);
541 
542  libMesh::out << "Refined mesh to "
543  << mesh.n_active_elem()
544  << " active elements and "
545  << equation_systems.n_active_dofs()
546  << " active dofs."
547  << std::endl;
548  }
549 
550  // Do one last solve if necessary
551  if (a_step == param.max_adaptivesteps)
552  {
553  system.solve();
554 
555  write_output(equation_systems, a_step, "primal", param);
556 
557  system.assemble_qoi_sides = true;
558 
559  SensitivityData sensitivities(qois, system, system.get_parameter_vector());
560 
561  // Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity
562  // function, so we have to set the adjoint_already_solved boolean to false
563  system.set_adjoint_already_solved(false);
564 
565  system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
566 
567  // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
568  system.set_adjoint_already_solved(true);
569 
570  GetPot infile_l_shaped("l-shaped.in");
571 
572  Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
573  Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
574  Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
575  Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
576 
577  libMesh::out << "Adaptive step "
578  << a_step
579  << ", we have "
580  << mesh.n_active_elem()
581  << " active elements and "
582  << equation_systems.n_active_dofs()
583  << " active dofs."
584  << std::endl;
585 
586  libMesh::out << "Sensitivity of QoI one to Parameter one is "
587  << sensitivity_QoI_0_0_computed
588  << std::endl;
589 
590  libMesh::out << "Sensitivity of QoI one to Parameter two is "
591  << sensitivity_QoI_0_1_computed
592  << std::endl;
593 
594  libMesh::out << "The error in sensitivity QoI_0_0 is "
595  << std::setprecision(17)
596  << std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact
597  << std::endl;
598 
599  libMesh::out << "The error in sensitivity QoI_0_1 is "
600  << std::setprecision(17)
601  << std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact
602  << std::endl
603  << std::endl;
604 
605  // Hard coded asserts to ensure that the actual numbers we are getting are what they should be
606  libmesh_assert_less(std::abs((sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact), 2.e-4);
607  libmesh_assert_less(std::abs((sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact), 2.e-4);
608 
609  // Let's do a forward sensitivity solve too, unless we're
610  // told to skip it for backwards compatibility with old
611  // performance benchmarks.
612  const bool forward_sensitivity = infile("--forward_sensitivity", true);
613 
614  // Don't confuse PETSc with our custom GetPot's arguments
616 
617  if (forward_sensitivity)
618  {
619  // This will require two linear solves (one per parameter)
620  // rather than the adjoint sensitivity's one, but it's useful
621  // for regression testing.
622  SensitivityData forward_sensitivities(qois, system, system.get_parameter_vector());
623  system.forward_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), forward_sensitivities);
624 
625  libmesh_assert_less(std::abs((forward_sensitivities[0][0] - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact), 2.e-4);
626  libmesh_assert_less(std::abs((forward_sensitivities[0][1] - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact), 2.e-4);
627 
628  // These should be the same linearization, just calculated
629  // different ways with different roundoff error
630  libmesh_assert_less
631  (std::abs((forward_sensitivities[0][0] - sensitivity_QoI_0_0_computed)/sensitivity_QoI_0_0_computed), TOLERANCE);
632  libmesh_assert_less
633  (std::abs((forward_sensitivities[0][1] - sensitivity_QoI_0_1_computed)/sensitivity_QoI_0_1_computed), TOLERANCE);
634 
635  libMesh::out << "The error in forward calculation of sensitivity QoI_0_0 is "
636  << std::setprecision(17)
637  << std::abs(forward_sensitivities[0][0] - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact
638  << std::endl;
639 
640  libMesh::out << "The error in forward calculation of sensitivity QoI_0_1 is "
641  << std::setprecision(17)
642  << std::abs(forward_sensitivities[0][1] - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact
643  << std::endl
644  << std::endl;
645 
646 
647  }
648 
649  NumericVector<Number> & primal_solution = *system.solution;
650  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
651  primal_solution.swap(dual_solution_0);
652  write_output(equation_systems, a_step, "adjoint_0", param);
653 
654  primal_solution.swap(dual_solution_0);
655  }
656  }
657 
658  libMesh::err << '[' << system.processor_id()
659  << "] Completing output."
660  << std::endl;
661 
662 #endif // #ifndef LIBMESH_ENABLE_AMR
663 
664  // All done.
665  return 0;
666 }
unsigned int nelem_target
Definition: femparameters.h:60
OStreamProxy err
bool print_solution_norms
This is the EquationSystems class.
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
virtual dof_id_type n_active_elem() const =0
Real verify_analytic_jacobians
If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calc...
Definition: fem_system.h:215
int main(int argc, char **argv)
Definition: adjoints_ex2.C:316
void set_adjoint_already_solved(bool setting)
Setter for the adjoint_already_solved boolean.
Definition: system.h:412
bool analytic_jacobians
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=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
std::unique_ptr< ErrorEstimator > build_error_estimator(FEMParameters &param)
Definition: adjoints_ex2.C:287
static constexpr Real TOLERANCE
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
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Definition: diff_system.h:359
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
void adjust_linear_solver(LinearSolver< Number > &linear_solver)
Definition: adjoints_ex2.C:169
std::string & fe_family()
Definition: L-shaped.h:24
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
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
MeshBase & mesh
This class provides an interface to Eigen iterative solvers that is compatible with the libMesh Linea...
bool print_residual_norms
libMesh::Real relative_step_tolerance
std::string indicator_type
bool & analytic_jacobians()
Definition: L-shaped.h:26
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
The libMesh namespace provides an interface to certain functionality in the library.
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
libMesh::Real refine_fraction
Definition: femparameters.h:62
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:364
virtual std::unique_ptr< DiffSolver > & diff_solver()
An implicit linear or nonlinear solver to use at each timestep.
Definition: time_solver.h:220
This is the MeshBase class.
Definition: mesh_base.h:74
std::vector< unsigned int > fe_order
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
virtual void forward_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) override
Solves for the derivative of each of the system&#39;s quantities of interest q in qoi[qoi_indices] with r...
SolverPackage default_solver_package()
Definition: libmesh.C:1050
libMesh::Real coarsen_threshold
Definition: femparameters.h:62
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 &...)
Implements (adaptive) mesh refinement algorithms for a MeshBase.
unsigned int max_linear_iterations
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
Data structure for holding completed parameter sensitivity calculations.
virtual void adjoint_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) override
Solves for the derivative of each of the system&#39;s quantities of interest q in qoi[qoi_indices] with r...
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Definition: mesh_base.C:1489
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
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
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.
double initial_linear_tolerance
libMesh::Real global_tolerance
Definition: femparameters.h:61
void attach_qoi(DifferentiableQoI *qoi_in)
Attach external QoI object.
Definition: diff_system.C:279
void set_solver_type(const SolverType st)
Sets the type of solver to use.
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
unsigned int & fe_order()
Definition: L-shaped.h:25
void add_command_line_names(const GetPot &getpot)
Merge a GetPot object&#39;s requested names into the set of queried command-line names.
Definition: libmesh.C:906
virtual LinearSolver< Number > * get_linear_solver() const override
Definition: diff_system.C:163
bool print_jacobian_norms
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
void adjust_linear_solvers(LaplaceSystem &system)
Definition: adjoints_ex2.C:184
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:99
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time, const std::set< std::string > *system_names=nullptr)
Writes out the solution at a specific timestep.
Definition: exodusII_io.C:1995
void write_output(EquationSystems &es, unsigned int a_step, std::string solution_type, FEMParameters &param)
Definition: adjoints_ex2.C:110
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
const MeshBase & get_mesh() const
libMesh::Real min_step_length
std::size_t n_active_dofs() const
void all_second_order(const bool full_ordered=true)
Calls the range-based version of this function with a range consisting of all elements in the mesh...
Definition: mesh_base.C:1535
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:1193
unsigned int max_nonlinear_iterations
double linear_tolerance_multiplier
virtual void solve() override
Invokes the solver associated with the system.
Definition: fem_system.C:1076
virtual void init()
Initialize all the systems.
std::unique_ptr< MeshRefinement > build_mesh_refinement(MeshBase &mesh, FEMParameters &param)
Definition: adjoints_ex2.C:264
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
processor_id_type processor_id() const
ParameterVector & get_parameter_vector()
Definition: L-shaped.h:37
void set_system_parameters(LaplaceSystem &system, FEMParameters &param)
Definition: adjoints_ex2.C:201
libMesh::Real coarsen_fraction
Definition: femparameters.h:62
TimeSolver & get_time_solver()
Definition: diff_system.h:454
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.