1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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
64 // General libMesh includes
65 #include "libmesh/equation_systems.h"
66 #include "libmesh/error_vector.h"
67 #include "libmesh/mesh.h"
68 #include "libmesh/mesh_refinement.h"
69 #include "libmesh/newton_solver.h"
70 #include "libmesh/numeric_vector.h"
71 #include "libmesh/steady_solver.h"
72 #include "libmesh/system_norm.h"
73 #include "libmesh/auto_ptr.h" // libmesh_make_unique
74 #include "libmesh/enum_solver_package.h"
75
76 // Sensitivity Calculation related includes
77 #include "libmesh/parameter_vector.h"
78 #include "libmesh/sensitivity_data.h"
79
80 // Error Estimator includes
81 #include "libmesh/kelly_error_estimator.h"
82 #include "libmesh/patch_recovery_error_estimator.h"
83
84 // Adjoint Related includes
85 #include "libmesh/adjoint_residual_error_estimator.h"
86 #include "libmesh/qoi_set.h"
87
88 // libMesh I/O includes
89 #include "libmesh/getpot.h"
90 #include "libmesh/gmv_io.h"
91 #include "libmesh/exodusII_io.h"
92
93 // Local includes
94 #include "femparameters.h"
95 #include "L-shaped.h"
96 #include "L-qoi.h"
97
98 // Bring in everything from the libMesh namespace
99 using namespace libMesh;
100
101 // Local function declarations
102
103 // Number output files, the files are give a prefix of primal or adjoint_i depending on
104 // whether the output is the primal solution or the dual solution for the ith QoI
105
106 // Write gmv output
write_output(EquationSystems & es,unsigned int a_step,std::string solution_type,FEMParameters & param)107 void write_output(EquationSystems & es,
108 unsigned int a_step, // The adaptive step count
109 std::string solution_type, // primal or adjoint solve
110 FEMParameters & param)
111 {
112 // Ignore parameters when there are no output formats available.
113 libmesh_ignore(es, a_step, solution_type, param);
114
115 #ifdef LIBMESH_HAVE_GMV
116 if (param.output_gmv)
117 {
118 MeshBase & mesh = es.get_mesh();
119
120 std::ostringstream file_name_gmv;
121 file_name_gmv << solution_type
122 << ".out.gmv."
123 << std::setw(2)
124 << std::setfill('0')
125 << std::right
126 << a_step;
127
128 GMVIO(mesh).write_equation_systems(file_name_gmv.str(), es);
129 }
130 #endif
131
132 #ifdef LIBMESH_HAVE_EXODUS_API
133 if (param.output_exodus)
134 {
135 MeshBase & mesh = es.get_mesh();
136
137 // We write out one file per adaptive step. The files are named in
138 // the following way:
139 // foo.e
140 // foo.e-s002
141 // foo.e-s003
142 // ...
143 // so that, if you open the first one with Paraview, it actually
144 // opens the entire sequence of adapted files.
145 std::ostringstream file_name_exodus;
146
147 file_name_exodus << solution_type << ".e";
148 if (a_step > 0)
149 file_name_exodus << "-s"
150 << std::setw(3)
151 << std::setfill('0')
152 << std::right
153 << a_step + 1;
154
155 // We write each adaptive step as a pseudo "time" step, where the
156 // time simply matches the (1-based) adaptive step we are on.
157 ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
158 es,
159 1,
160 /*time=*/a_step + 1);
161 }
162 #endif
163 }
164
165 // Set the parameters for the nonlinear and linear solvers to be used during the simulation
set_system_parameters(LaplaceSystem & system,FEMParameters & param)166 void set_system_parameters(LaplaceSystem & system, FEMParameters & param)
167 {
168 // Use analytical jacobians?
169 system.analytic_jacobians() = param.analytic_jacobians;
170
171 // Verify analytic jacobians against numerical ones?
172 system.verify_analytic_jacobians = param.verify_analytic_jacobians;
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 // More desperate debugging options
179 system.print_solution_norms = param.print_solution_norms;
180 system.print_solutions = param.print_solutions;
181 system.print_residual_norms = param.print_residual_norms;
182 system.print_residuals = param.print_residuals;
183 system.print_jacobian_norms = param.print_jacobian_norms;
184 system.print_jacobians = param.print_jacobians;
185
186 // No transient time solver
187 system.time_solver = libmesh_make_unique<SteadySolver>(system);
188
189 // Nonlinear solver options
190 {
191 NewtonSolver * solver = new NewtonSolver(system);
192 system.time_solver->diff_solver() = std::unique_ptr<DiffSolver>(solver);
193
194 solver->quiet = param.solver_quiet;
195 solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
196 solver->minsteplength = param.min_step_length;
197 solver->relative_step_tolerance = param.relative_step_tolerance;
198 solver->relative_residual_tolerance = param.relative_residual_tolerance;
199 solver->require_residual_reduction = param.require_residual_reduction;
200 solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
201 if (system.time_solver->reduce_deltat_on_diffsolver_failure)
202 {
203 solver->continue_after_max_iterations = true;
204 solver->continue_after_backtrack_failure = true;
205 }
206
207 // And the linear solver options
208 solver->max_linear_iterations = param.max_linear_iterations;
209 solver->initial_linear_tolerance = param.initial_linear_tolerance;
210 solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
211 }
212 }
213
214 // Build the mesh refinement object and set parameters for refining/coarsening etc
215
216 #ifdef LIBMESH_ENABLE_AMR
217
build_mesh_refinement(MeshBase & mesh,FEMParameters & param)218 std::unique_ptr<MeshRefinement> build_mesh_refinement(MeshBase & mesh,
219 FEMParameters & param)
220 {
221 MeshRefinement * mesh_refinement = new MeshRefinement(mesh);
222 mesh_refinement->coarsen_by_parents() = true;
223 mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
224 mesh_refinement->nelem_target() = param.nelem_target;
225 mesh_refinement->refine_fraction() = param.refine_fraction;
226 mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
227 mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
228
229 return std::unique_ptr<MeshRefinement>(mesh_refinement);
230 }
231
232 #endif // LIBMESH_ENABLE_AMR
233
234 // This is where we declare the error estimators to be built and used for
235 // mesh refinement. The adjoint residual estimator needs two estimators.
236 // One for the forward component of the estimate and one for the adjoint
237 // weighting factor. Here we use the Patch Recovery indicator to estimate both the
238 // forward and adjoint weights. The H1 seminorm component of the error is used
239 // as dictated by the weak form the Laplace equation.
240
build_error_estimator(FEMParameters & param)241 std::unique_ptr<ErrorEstimator> build_error_estimator(FEMParameters & param)
242 {
243 if (param.indicator_type == "kelly")
244 {
245 libMesh::out << "Using Kelly Error Estimator" << std::endl;
246
247 return libmesh_make_unique<KellyErrorEstimator>();
248 }
249 else if (param.indicator_type == "adjoint_residual")
250 {
251 libMesh::out << "Using Adjoint Residual Error Estimator with Patch Recovery Weights" << std::endl << std::endl;
252
253 AdjointResidualErrorEstimator * adjoint_residual_estimator = new AdjointResidualErrorEstimator;
254
255 adjoint_residual_estimator->error_plot_suffix = "error.gmv";
256
257 PatchRecoveryErrorEstimator * p1 = new PatchRecoveryErrorEstimator;
258 adjoint_residual_estimator->primal_error_estimator().reset(p1);
259
260 PatchRecoveryErrorEstimator * p2 = new PatchRecoveryErrorEstimator;
261 adjoint_residual_estimator->dual_error_estimator().reset(p2);
262
263 adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
264
265 adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
266
267 return std::unique_ptr<ErrorEstimator>(adjoint_residual_estimator);
268 }
269 else
270 libmesh_error_msg("Unknown indicator_type = " << param.indicator_type);
271 }
272
273 // The main program.
main(int argc,char ** argv)274 int main (int argc, char ** argv)
275 {
276 // Initialize libMesh.
277 LibMeshInit init (argc, argv);
278
279 // This example requires a linear solver package.
280 libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
281 "--enable-petsc, --enable-trilinos, or --enable-eigen");
282
283 // Skip adaptive examples on a non-adaptive libMesh build
284 #ifndef LIBMESH_ENABLE_AMR
285 libmesh_example_requires(false, "--enable-amr");
286 #else
287
288 libMesh::out << "Started " << argv[0] << std::endl;
289
290 // Make sure the general input file exists, and parse it
291 {
292 std::ifstream i("general.in");
293 libmesh_error_msg_if(!i, '[' << init.comm().rank() << "] Can't find general.in; exiting early.");
294 }
295 GetPot infile("general.in");
296
297 // Read in parameters from the input file
298 FEMParameters param(init.comm());
299 param.read(infile);
300
301 // Skip this default-2D example if libMesh was compiled as 1D-only.
302 libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
303
304 // Create a mesh, with dimension to be overridden later, distributed
305 // across the default MPI communicator.
306 Mesh mesh(init.comm());
307
308 // And an object to refine it
309 std::unique_ptr<MeshRefinement> mesh_refinement =
310 build_mesh_refinement(mesh, param);
311
312 // And an EquationSystems to run on it
313 EquationSystems equation_systems (mesh);
314
315 libMesh::out << "Reading in and building the mesh" << std::endl;
316
317 // Read in the mesh
318 mesh.read(param.domainfile.c_str());
319 // Make all the elements of the mesh second order so we can compute
320 // with a higher order basis
321 mesh.all_second_order();
322
323 // Create a mesh refinement object to do the initial uniform refinements
324 // on the coarse grid read in from lshaped.xda
325 MeshRefinement initial_uniform_refinements(mesh);
326 initial_uniform_refinements.uniformly_refine(param.coarserefinements);
327
328 libMesh::out << "Building system" << std::endl;
329
330 // Build the FEMSystem
331 LaplaceSystem & system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
332
333 QoISet qois;
334
335 std::vector<unsigned int> qoi_indices;
336 qoi_indices.push_back(0);
337 qois.add_indices(qoi_indices);
338
339 qois.set_weight(0, 0.5);
340
341 // Put some scope here to test that the cloning is working right
342 {
343 LaplaceQoI qoi;
344 system.attach_qoi(&qoi);
345 }
346
347 // Set its parameters
348 set_system_parameters(system, param);
349
350 libMesh::out << "Initializing systems" << std::endl;
351
352 equation_systems.init ();
353
354 // Print information about the mesh and system to the screen.
355 mesh.print_info();
356 equation_systems.print_info();
357
358 {
359 // Adaptively solve the timestep
360 unsigned int a_step = 0;
361 for (; a_step != param.max_adaptivesteps; ++a_step)
362 {
363 // We can't adapt to both a tolerance and a
364 // target mesh size
365 if (param.global_tolerance != 0.)
366 libmesh_assert_equal_to (param.nelem_target, 0);
367 // If we aren't adapting to a tolerance we need a
368 // target mesh size
369 else
370 libmesh_assert_greater (param.nelem_target, 0);
371
372 // Solve the forward problem
373 system.solve();
374
375 // Write out the computed primal solution
376 write_output(equation_systems, a_step, "primal", param);
377
378 // Get a pointer to the primal solution vector
379 NumericVector<Number> & primal_solution = *system.solution;
380
381 // A SensitivityData object to hold the qois and parameters
382 SensitivityData sensitivities(qois, system, system.get_parameter_vector());
383
384 // Make sure we get the contributions to the adjoint RHS from the sides
385 system.assemble_qoi_sides = true;
386
387 // Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity
388 // function, so we have to set the adjoint_already_solved boolean to false
389 system.set_adjoint_already_solved(false);
390
391 // Compute the sensitivities
392 system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
393
394 // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
395 system.set_adjoint_already_solved(true);
396
397 GetPot infile_l_shaped("l-shaped.in");
398
399 Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
400 Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
401 Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
402 Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
403
404 libMesh::out << "Adaptive step "
405 << a_step
406 << ", we have "
407 << mesh.n_active_elem()
408 << " active elements and "
409 << equation_systems.n_active_dofs()
410 << " active dofs."
411 << std::endl;
412
413 libMesh::out << "Sensitivity of QoI one to Parameter one is "
414 << sensitivity_QoI_0_0_computed
415 << std::endl;
416 libMesh::out << "Sensitivity of QoI one to Parameter two is "
417 << sensitivity_QoI_0_1_computed
418 << std::endl;
419
420 libMesh::out << "The relative error in sensitivity QoI_0_0 is "
421 << std::setprecision(17)
422 << std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact) / std::abs(sensitivity_QoI_0_0_exact)
423 << std::endl;
424
425 libMesh::out << "The relative error in sensitivity QoI_0_1 is "
426 << std::setprecision(17)
427 << std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact) / std::abs(sensitivity_QoI_0_1_exact)
428 << std::endl
429 << std::endl;
430
431 // Get a pointer to the solution vector of the adjoint problem for QoI 0
432 NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
433
434 // Swap the primal and dual solutions so we can write out the adjoint solution
435 primal_solution.swap(dual_solution_0);
436 write_output(equation_systems, a_step, "adjoint_0", param);
437
438 // Swap back
439 primal_solution.swap(dual_solution_0);
440
441 // We have to refine either based on reaching an error tolerance or
442 // a number of elements target, which should be verified above
443 // Otherwise we flag elements by error tolerance or nelem target
444
445 // Uniform refinement
446 if (param.refine_uniformly)
447 {
448 libMesh::out << "Refining Uniformly" << std::endl << std::endl;
449
450 mesh_refinement->uniformly_refine(1);
451 }
452 // Adaptively refine based on reaching an error tolerance
453 else if (param.global_tolerance >= 0. && param.nelem_target == 0.)
454 {
455 // Now we construct the data structures for the mesh refinement process
456 ErrorVector error;
457
458 // Build an error estimator object
459 std::unique_ptr<ErrorEstimator> error_estimator =
460 build_error_estimator(param);
461
462 // Estimate the error in each element using the Adjoint Residual or Kelly error estimator
463 error_estimator->estimate_error(system, error);
464
465 mesh_refinement->flag_elements_by_error_tolerance (error);
466
467 mesh_refinement->refine_and_coarsen_elements();
468 }
469 // Adaptively refine based on reaching a target number of elements
470 else
471 {
472 // Now we construct the data structures for the mesh refinement process
473 ErrorVector error;
474
475 // Build an error estimator object
476 std::unique_ptr<ErrorEstimator> error_estimator =
477 build_error_estimator(param);
478
479 // Estimate the error in each element using the Adjoint Residual or Kelly error estimator
480 error_estimator->estimate_error(system, error);
481
482 if (mesh.n_active_elem() >= param.nelem_target)
483 {
484 libMesh::out<<"We reached the target number of elements."<<std::endl <<std::endl;
485 break;
486 }
487
488 mesh_refinement->flag_elements_by_nelem_target (error);
489
490 mesh_refinement->refine_and_coarsen_elements();
491 }
492
493 // Dont forget to reinit the system after each adaptive refinement !
494 equation_systems.reinit();
495
496 libMesh::out << "Refined mesh to "
497 << mesh.n_active_elem()
498 << " active elements and "
499 << equation_systems.n_active_dofs()
500 << " active dofs."
501 << std::endl;
502 }
503
504 // Do one last solve if necessary
505 if (a_step == param.max_adaptivesteps)
506 {
507 system.solve();
508
509 write_output(equation_systems, a_step, "primal", param);
510
511 NumericVector<Number> & primal_solution = *system.solution;
512
513 SensitivityData sensitivities(qois, system, system.get_parameter_vector());
514
515 system.assemble_qoi_sides = true;
516
517 // Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity
518 // function, so we have to set the adjoint_already_solved boolean to false
519 system.set_adjoint_already_solved(false);
520
521 system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
522
523 // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
524 system.set_adjoint_already_solved(true);
525
526 GetPot infile_l_shaped("l-shaped.in");
527
528 Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
529 Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
530 Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
531 Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
532
533 libMesh::out << "Adaptive step "
534 << a_step
535 << ", we have "
536 << mesh.n_active_elem()
537 << " active elements and "
538 << equation_systems.n_active_dofs()
539 << " active dofs."
540 << std::endl;
541
542 libMesh::out << "Sensitivity of QoI one to Parameter one is "
543 << sensitivity_QoI_0_0_computed
544 << std::endl;
545
546 libMesh::out << "Sensitivity of QoI one to Parameter two is "
547 << sensitivity_QoI_0_1_computed
548 << std::endl;
549
550 libMesh::out << "The error in sensitivity QoI_0_0 is "
551 << std::setprecision(17)
552 << std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact
553 << std::endl;
554
555 libMesh::out << "The error in sensitivity QoI_0_1 is "
556 << std::setprecision(17)
557 << std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact
558 << std::endl
559 << std::endl;
560
561 // Hard coded asserts to ensure that the actual numbers we are getting are what they should be
562 libmesh_assert_less(std::abs((sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact), 2.e-4);
563 libmesh_assert_less(std::abs((sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact), 2.e-4);
564
565 NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
566
567 primal_solution.swap(dual_solution_0);
568 write_output(equation_systems, a_step, "adjoint_0", param);
569
570 primal_solution.swap(dual_solution_0);
571 }
572 }
573
574 libMesh::err << '[' << system.processor_id()
575 << "] Completing output."
576 << std::endl;
577
578 #endif // #ifndef LIBMESH_ENABLE_AMR
579
580 // All done.
581 return 0;
582 }
583