libMesh
adjoints_ex3.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // <h1>Adjoints Example 3 - Stokes flow and Convection-Diffusion</h1>
20 // \author Vikram Garg
21 // \date 2012
22 //
23 // We solve a coupled Stokes + Convection Diffusion system in an H channel geometry
24 // with 2 inlets and 2 outlets. The QoI is the species
25 // flux from left outlet
26 
27 // Channel Geometry:
28 // Wall
29 // ----------------------------------------------------------------
30 // 0 1
31 // ---------------------------- -------------------------------
32 // | |
33 // Wall | | Wall
34 // | |
35 // ----------------------------- -------------------------------
36 // 2 2
37 // -----------------------------------------------------------------
38 // Wall
39 
40 // The equations governing this flow are:
41 // Stokes: -VectorLaplacian(velocity) + grad(pressure) = vector(0)
42 // Convection-Diffusion: - dot(velocity, grad(concentration) ) + Laplacian(concentration) = 0
43 
44 // The boundary conditions are:
45 // u_1(0) = -(y-2)*(y-3), u_2(0) = 0 ; u_1(1) = (y-2)*(y-3), u_2(1) = 0 ;
46 // u_1(walls) = 0, u_2(walls) = 0;
47 // C(0) = 1 ; C(1) = 0;
48 // grad(C) dot n (walls) = 0 ;
49 // grad(C) dot n (2) = 0 ; grad(C) dot n (3) = 0
50 
51 // The QoI is:
52 // Q((u,v), C) = integral_{left_outlet} - u * C ds
53 
54 // The complete equal order adjoint QoI error estimate is: (Notation for derivatives: grad(C) = C,1 + C,2)
55 // Q(u) - Q(u_h) \leq
56 // |e(u_1)|_{H1} |e(u_1^*)|_{H1} + |e(u_2)|_{H1} |e(u_2^*)|_{H1} + (1/Pe) * |e(C)|_{H1} |e(C^*)|_{H1} +
57 // ||e(u_1,1)||_{L2} ||e(p^*)||_{L2} + ||e(u_2,2)||_{L2} ||e(p^*)||_{L2} + ||e(u_1,1^*)||_{L2} ||e(p)||_{L2} + ||e(u_2,2^*)||_{L2} ||e(p)||_{L2} +
58 // ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} + ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2}
59 // ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} + ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2}
60 // = error_non_pressure + error_with_pressure + error_convection_diffusion_x + error_convection_diffusion_y
61 
62 // C++ includes
63 #include <iostream>
64 #include <sys/time.h>
65 #include <iomanip>
66 
67 // Libmesh includes
68 #include "libmesh/equation_systems.h"
69 
70 #include "libmesh/twostep_time_solver.h"
71 #include "libmesh/euler_solver.h"
72 #include "libmesh/euler2_solver.h"
73 #include "libmesh/steady_solver.h"
74 
75 #include "libmesh/newton_solver.h"
76 #include "libmesh/numeric_vector.h"
77 #include "libmesh/petsc_diff_solver.h"
78 
79 #include "libmesh/mesh.h"
80 #include "libmesh/mesh_tools.h"
81 #include "libmesh/mesh_base.h"
82 #include "libmesh/mesh_refinement.h"
83 
84 #include "libmesh/system.h"
85 #include "libmesh/system_norm.h"
86 
87 #include "libmesh/adjoint_residual_error_estimator.h"
88 #include "libmesh/const_fem_function.h"
89 #include "libmesh/error_vector.h"
90 #include "libmesh/fem_function_base.h"
91 #include "libmesh/getpot.h"
92 #include "libmesh/gmv_io.h"
93 #include "libmesh/exodusII_io.h"
94 #include "libmesh/kelly_error_estimator.h"
95 #include "libmesh/parameter_vector.h"
96 #include "libmesh/patch_recovery_error_estimator.h"
97 #include "libmesh/petsc_vector.h"
98 #include "libmesh/sensitivity_data.h"
99 #include "libmesh/tecplot_io.h"
100 #include "libmesh/uniform_refinement_estimator.h"
101 #include "libmesh/qoi_set.h"
102 #include "libmesh/weighted_patch_recovery_error_estimator.h"
103 
104 // Local includes
105 #include "coupled_system.h"
106 #include "domain.h"
107 #include "initial.h"
108 #include "femparameters.h"
109 #include "mysystems.h"
110 #include "output.h"
111 #include "H-qoi.h"
112 
113 
114 // Bring in everything from the libMesh namespace
115 using namespace libMesh;
116 
117 // Number output files
118 
119 std::string numbered_filename(unsigned int t_step, // The timestep count
120  unsigned int a_step, // The adaptive step count
121  std::string solution_type, // primal or adjoint solve
122  std::string type,
123  std::string extension,
124  FEMParameters & param)
125 {
126  std::ostringstream file_name;
127  file_name << solution_type
128  << ".out."
129  << type
130  << '.'
131  << std::setw(3)
132  << std::setfill('0')
133  << std::right
134  << t_step
135  << '.'
136  << std::setw(2)
137  << std::setfill('0')
138  << std::right
139  << a_step
140  << '.'
141  << extension;
142 
143  if (param.output_bz2)
144  file_name << ".bz2";
145  else if (param.output_gz)
146  file_name << ".gz";
147  return file_name.str();
148 }
149 
150 // Possibly write tecplot, gmv, xda/xdr, and exodus output files.
152  unsigned int t_step, // The timestep count
153  unsigned int a_step, // The adaptive step count
154  std::string solution_type, // primal or adjoint solve
155  FEMParameters & param)
156 {
157  MeshBase & mesh = es.get_mesh();
158 
159 #ifdef LIBMESH_HAVE_GMV
160  if (param.output_gmv)
161  {
162  std::ostringstream file_name_gmv;
163  file_name_gmv << solution_type
164  << ".out.gmv."
165  << std::setw(3)
166  << std::setfill('0')
167  << std::right
168  << t_step
169  << '.'
170  << std::setw(2)
171  << std::setfill('0')
172  << std::right
173  << a_step;
174 
176  (file_name_gmv.str(), es);
177  }
178 #endif
179 
180 #ifdef LIBMESH_HAVE_TECPLOT_API
181  if (param.output_tecplot)
182  {
183  std::ostringstream file_name_tecplot;
184  file_name_tecplot << solution_type
185  << ".out."
186  << std::setw(3)
187  << std::setfill('0')
188  << std::right
189  << t_step
190  << '.'
191  << std::setw(2)
192  << std::setfill('0')
193  << std::right
194  << a_step
195  << ".plt";
196 
198  (file_name_tecplot.str(), es);
199  }
200 #endif
201 
202  if (param.output_xda || param.output_xdr)
204  if (param.output_xda)
205  {
206  mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param));
207  es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xda", param),
210  }
211  if (param.output_xdr)
212  {
213  mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param));
214  es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param),
217  }
218 
219 #ifdef LIBMESH_HAVE_EXODUS_API
220  if (param.output_exodus)
221  {
222  // We write out one file per adaptive step. The files are named in
223  // the following way:
224  // foo.e
225  // foo.e-s002
226  // foo.e-s003
227  // ...
228  // so that, if you open the first one with Paraview, it actually
229  // opens the entire sequence of adapted files.
230  std::ostringstream file_name_exodus;
231 
232  file_name_exodus << solution_type << ".e";
233  if (a_step > 0)
234  file_name_exodus << "-s"
235  << std::setw(3)
236  << std::setfill('0')
237  << std::right
238  << a_step + 1;
239 
240  // We write each adaptive step as a pseudo "time" step, where the
241  // time simply matches the (1-based) adaptive step we are on.
242  ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
243  es,
244  1,
245  /*time=*/a_step + 1);
246  }
247 #endif
248 }
249 
251 {
252  // Only one processor needs to take care of headers.
253  if (libMesh::global_processor_id() != 0)
254  return;
255 
256  start_output(param.initial_timestep, "out_clocktime.m", "vector_clocktime");
257 
258  if (param.run_simulation)
259  {
260  start_output(param.initial_timestep, "out_activemesh.m", "table_activemesh");
261  start_output(param.initial_timestep, "out_solvesteps.m", "table_solvesteps");
262 
263  if (param.timesolver_tolerance)
264  {
265  start_output(param.initial_timestep, "out_time.m", "vector_time");
266  start_output(param.initial_timestep, "out_timesteps.m", "vector_timesteps");
267  }
268  if (param.steadystate_tolerance)
269  start_output(param.initial_timestep, "out_changerate.m", "vector_changerate");
270  }
271 }
272 
274  unsigned int a_step,
275  unsigned int newton_steps,
276  unsigned int krylov_steps,
277  unsigned int tv_sec,
278  unsigned int tv_usec)
279 {
280  MeshBase & mesh = es.get_mesh();
281  unsigned int n_active_elem = mesh.n_active_elem();
282  unsigned int n_active_dofs = es.n_active_dofs();
283 
284  if (mesh.processor_id() == 0)
285  {
286  // Write out the number of elements/dofs used
287  std::ofstream activemesh ("out_activemesh.m",
288  std::ios_base::app | std::ios_base::out);
289  activemesh.precision(17);
290  activemesh << (a_step + 1) << ' '
291  << n_active_elem << ' '
292  << n_active_dofs << std::endl;
293 
294  // Write out the number of solver steps used
295  std::ofstream solvesteps ("out_solvesteps.m",
296  std::ios_base::app | std::ios_base::out);
297  solvesteps.precision(17);
298  solvesteps << newton_steps << ' '
299  << krylov_steps << std::endl;
300 
301  // Write out the clock time
302  std::ofstream clocktime ("out_clocktime.m",
303  std::ios_base::app | std::ios_base::out);
304  clocktime.precision(17);
305  clocktime << tv_sec << '.' << tv_usec << std::endl;
306  }
307 }
308 
310 {
311  // Write footers on output .m files
312  if (libMesh::global_processor_id() == 0)
313  {
314  std::ofstream clocktime ("out_clocktime.m",
315  std::ios_base::app | std::ios_base::out);
316  clocktime << "];" << std::endl;
317 
318  if (param.run_simulation)
319  {
320  std::ofstream activemesh ("out_activemesh.m",
321  std::ios_base::app | std::ios_base::out);
322  activemesh << "];" << std::endl;
323 
324  std::ofstream solvesteps ("out_solvesteps.m",
325  std::ios_base::app | std::ios_base::out);
326  solvesteps << "];" << std::endl;
327 
328  if (param.timesolver_tolerance)
329  {
330  std::ofstream times ("out_time.m",
331  std::ios_base::app | std::ios_base::out);
332  times << "];" << std::endl;
333  std::ofstream timesteps ("out_timesteps.m",
334  std::ios_base::app | std::ios_base::out);
335  timesteps << "];" << std::endl;
336  }
337  if (param.steadystate_tolerance)
338  {
339  std::ofstream changerate ("out_changerate.m",
340  std::ios_base::app | std::ios_base::out);
341  changerate << "];" << std::endl;
342  }
343  }
344 
345  // We're done, so write out a file (for e.g. Make to check)
346  std::ofstream complete ("complete");
347  complete << "complete" << std::endl;
348  }
349 }
350 
351 #if defined(LIBMESH_HAVE_GMV) || defined(LIBMESH_HAVE_TECPLOT_API)
352 void write_error(EquationSystems & es,
353  ErrorVector & error,
354  unsigned int t_number,
355  unsigned int a_number,
356  FEMParameters & param,
357  std::string error_type)
358 #else
361  unsigned int,
362  unsigned int,
363  FEMParameters &,
364  std::string)
365 #endif
366 {
367 #ifdef LIBMESH_HAVE_GMV
368  if (param.write_gmv_error)
369  {
370  std::ostringstream error_gmv;
371  error_gmv << "error.gmv."
372  << std::setw(3)
373  << std::setfill('0')
374  << std::right
375  << a_number
376  << "."
377  << std::setw(2)
378  << std::setfill('0')
379  << std::right
380  << t_number
381  << error_type;
382 
383  error.plot_error(error_gmv.str(), es.get_mesh());
384  }
385 #endif
386 
387 #ifdef LIBMESH_HAVE_TECPLOT_API
388  if (param.write_tecplot_error)
389  {
390  std::ostringstream error_tecplot;
391  error_tecplot << "error.plt."
392  << std::setw(3)
393  << std::setfill('0')
394  << std::right
395  << a_number
396  << "."
397  << std::setw(2)
398  << std::setfill('0')
399  << std::right
400  << t_number
401  << error_type;
402 
403  error.plot_error(error_tecplot.str(), es.get_mesh());
404  }
405 #endif
406 }
407 
408 void read_output(EquationSystems & es,
409  unsigned int t_step,
410  unsigned int a_step,
411  std::string solution_type,
412  FEMParameters & param)
413 {
414  MeshBase & mesh = es.get_mesh();
415 
416  std::string file_name_mesh, file_name_soln;
417  // Look for ASCII files first
418  if (param.output_xda)
419  {
420  file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param);
421  file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xda", param);
422  }
423  else if (param.output_xdr)
424  {
425  file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param);
426  file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param);
427  }
428 
429  // Read in the mesh
430  mesh.read(file_name_mesh);
431 
432  // And the stored solution
433  es.read(file_name_soln, READ,
437 
438  // Put systems in a consistent state
439  for (unsigned int i = 0; i != es.n_systems(); ++i)
440  es.get_system<FEMSystem>(i).update();
441 
442  // Figure out the current time
443  Real current_time = 0., current_timestep = 0.;
444 
445  if (param.timesolver_tolerance)
446  {
447  std::ifstream times ("out_time.m");
448  std::ifstream timesteps ("out_timesteps.m");
449  if (times.is_open() && timesteps.is_open())
450  {
451  // Read headers
452  const unsigned int headersize = 25;
453  char header[headersize];
454  timesteps.getline (header, headersize);
455  if (strcmp(header, "vector_timesteps = [") != 0)
456  libmesh_error_msg("Bad header in out_timesteps.m:\n" << header);
457 
458  times.getline (header, headersize);
459  if (strcmp(header, "vector_time = [") != 0)
460  libmesh_error_msg("Bad header in out_time.m:\n" << header);
461 
462  // Read each timestep
463  for (unsigned int i = 0; i != t_step; ++i)
464  {
465  if (!times.good())
466  libmesh_error_msg("Error: File out_time.m is in non-good state.");
467  times >> current_time;
468  timesteps >> current_timestep;
469  }
470  // Remember to increment the last timestep; out_times.m
471  // lists each *start* time
472  current_time += current_timestep;
473  }
474  else
475  libmesh_error_msg("Error opening out_time.m or out_timesteps.m");
476  }
477  else
478  current_time = t_step * param.deltat;
479 
480  for (unsigned int i = 0; i != es.n_systems(); ++i)
481  es.get_system<FEMSystem>(i).time = current_time;
482 }
483 
484 void set_system_parameters(FEMSystem & system,
485  FEMParameters & param)
486 {
487  // Verify analytic jacobians against numerical ones?
489 
490  // More desperate debugging options
492  system.print_solutions = param.print_solutions;
494  system.print_residuals = param.print_residuals;
496  system.print_jacobians = param.print_jacobians;
497 
500 
501  // Solve this as a time-dependent or steady system
502  if (param.transient)
503  {
504  UnsteadySolver *innersolver;
505  if (param.timesolver_core == "euler2")
506  {
507  Euler2Solver *euler2solver =
508  new Euler2Solver(system);
509 
510  euler2solver->theta = param.timesolver_theta;
511  innersolver = euler2solver;
512  }
513  else if (param.timesolver_core == "euler")
514  {
515  EulerSolver *eulersolver =
516  new EulerSolver(system);
517 
518  eulersolver->theta = param.timesolver_theta;
519  innersolver = eulersolver;
520  }
521  else
522  libmesh_error_msg("Don't recognize core TimeSolver type: " << param.timesolver_core);
523 
524  if (param.timesolver_tolerance)
525  {
526  TwostepTimeSolver *timesolver =
527  new TwostepTimeSolver(system);
528 
529  timesolver->max_growth = param.timesolver_maxgrowth;
530  timesolver->target_tolerance = param.timesolver_tolerance;
531  timesolver->upper_tolerance = param.timesolver_upper_tolerance;
532  timesolver->component_norm = SystemNorm(param.timesolver_norm);
533 
534  timesolver->core_time_solver.reset(innersolver);
535  system.time_solver.reset(timesolver);
536  }
537  else
538  system.time_solver.reset(innersolver);
539  }
540  else
541  system.time_solver.reset(new SteadySolver(system));
542 
543  system.time_solver->reduce_deltat_on_diffsolver_failure =
544  param.deltat_reductions;
545  system.time_solver->quiet = param.time_solver_quiet;
546 
547  // Set the time stepping options
548  system.deltat = param.deltat;
549 
550  // And the integration options
552 
553  // And the nonlinear solver options
554  if (param.use_petsc_snes)
555  {
556 #ifdef LIBMESH_HAVE_PETSC
557  PetscDiffSolver *solver = new PetscDiffSolver(system);
558  system.time_solver->diff_solver().reset(solver);
559 #else
560  libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
561 #endif
562  }
563  else
564  {
565  NewtonSolver *solver = new NewtonSolver(system);
566  system.time_solver->diff_solver().reset(solver);
567 
568  solver->quiet = param.solver_quiet;
569  solver->verbose = !param.solver_quiet;
571  solver->minsteplength = param.min_step_length;
577  if (system.time_solver->reduce_deltat_on_diffsolver_failure)
578  {
579  solver->continue_after_max_iterations = true;
580  solver->continue_after_backtrack_failure = true;
581  }
582 
583  // And the linear solver options
587  }
588 }
589 
590 #ifdef LIBMESH_ENABLE_AMR
591 
593  FEMParameters & param)
594 {
595  MeshRefinement * mesh_refinement = new MeshRefinement(mesh);
596  mesh_refinement->coarsen_by_parents() = true;
597  mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
598  mesh_refinement->nelem_target() = param.nelem_target;
599  mesh_refinement->refine_fraction() = param.refine_fraction;
600  mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
601  mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
602 
603  return UniquePtr<MeshRefinement>(mesh_refinement);
604 }
605 
606 #endif // LIBMESH_ENABLE_AMR
607 
608 // This function builds the Kelly error indicator. This indicator can be used
609 // for comparisons of adjoint and non-adjoint based error indicators
611 {
613 }
614 
615 // Functions to build the adjoint based error indicators
616 // The error_non_pressure and error_pressure contributions are estimated using
617 // the build_error_estimator_component_wise function below
619 build_error_estimator_component_wise (FEMParameters & param,
620  std::vector<std::vector<Real>> & term_weights,
621  std::vector<FEMNormType> & primal_error_norm_type,
622  std::vector<FEMNormType> & dual_error_norm_type)
623 {
624  AdjointResidualErrorEstimator * adjoint_residual_estimator = new AdjointResidualErrorEstimator;
625 
626  // Both the primal and dual weights are going to be estimated using the patch recovery error estimator
628  adjoint_residual_estimator->primal_error_estimator().reset(p1);
629 
631  adjoint_residual_estimator->dual_error_estimator().reset(p2);
632 
633  // Set the boolean for specifying whether we are reusing patches while building the patch recovery estimates
634  p1->set_patch_reuse(param.patch_reuse);
635  p2->set_patch_reuse(param.patch_reuse);
636 
637  // Using the user filled error norm type vector, we pass the type of norm to be used for
638  // the error in each variable, we can have different types of norms for the primal and
639  // dual variables
640  std::size_t size = primal_error_norm_type.size();
641 
642  libmesh_assert_equal_to (size, dual_error_norm_type.size());
643  for (std::size_t i = 0; i != size; ++i)
644  {
645  adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
646  adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
647  }
648 
649  // Now we set the right weights for each term in the error estimate, using the user provided
650  // term_weights matrix
651  libmesh_assert_equal_to (size, term_weights.size());
652  for (std::size_t i = 0; i != term_weights.size(); ++i)
653  {
654  libmesh_assert_equal_to (size, term_weights[i].size());
655  adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
656  for (std::size_t j = 0; j != size; ++j)
657  if (i != j)
658  adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
659  }
660 
661  return UniquePtr<ErrorEstimator>(adjoint_residual_estimator);
662 }
663 
664 // The error_convection_diffusion_x and error_convection_diffusion_y are the nonlinear contributions which
665 // are computed using the build_weighted_error_estimator_component_wise below
667 build_weighted_error_estimator_component_wise (FEMParameters & param,
668  std::vector<std::vector<Real>> & term_weights,
669  std::vector<FEMNormType> & primal_error_norm_type,
670  std::vector<FEMNormType> & dual_error_norm_type,
671  std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions)
672 {
673  AdjointResidualErrorEstimator * adjoint_residual_estimator = new AdjointResidualErrorEstimator;
674 
675  // Using the user filled error norm type vector, we pass the type of norm to be used for
676  // the error in each variable, we can have different types of norms for the primal and
677  // dual variables
678 
680  adjoint_residual_estimator->primal_error_estimator().reset(p1);
681 
683  adjoint_residual_estimator->dual_error_estimator().reset(p2);
684 
685  p1->set_patch_reuse(param.patch_reuse);
686  p2->set_patch_reuse(param.patch_reuse);
687 
688  // This is the critical difference with the build_error_estimate_component_wise
689  p1->weight_functions.clear();
690 
691  // We pass the pointers to the user specified weight functions to the patch recovery
692  // error estimator objects declared above
693  std::size_t size = primal_error_norm_type.size();
694 
695  libmesh_assert(coupled_system_weight_functions.size() == size);
696  libmesh_assert(dual_error_norm_type.size() == size);
697 
698  for (unsigned int i = 0; i != size; ++i)
699  {
700  p1->weight_functions.push_back(coupled_system_weight_functions[i]);
701  adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
702  adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
703  }
704 
705  // Now we set the right weights for each term in the error estimate, using the user provided
706  // term_weights matrix
707  libmesh_assert_equal_to (size, term_weights.size());
708  for (std::size_t i = 0; i != term_weights.size(); ++i)
709  {
710  libmesh_assert_equal_to (size, term_weights[i].size());
711  adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
712  for (std::size_t j = 0; j != size; ++j)
713  if (i != j)
714  adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
715  }
716 
717  return UniquePtr<ErrorEstimator>(adjoint_residual_estimator);
718 }
719 
720 // The main program.
721 int main (int argc, char ** argv)
722 {
723  // Initialize libMesh.
724  LibMeshInit init (argc, argv);
725 
726  // Skip adaptive examples on a non-adaptive libMesh build
727 #ifndef LIBMESH_ENABLE_AMR
728  libmesh_example_requires(false, "--enable-amr");
729 #else
730 
731  // This doesn't converge with Eigen BICGSTAB or with Trilinos for some reason...
732  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
733 
734  libMesh::out << "Started " << argv[0] << std::endl;
735 
736  // Make sure the general input file exists, and parse it
737  {
738  std::ifstream i("general.in");
739  if (!i)
740  libmesh_error_msg('[' << init.comm().rank() << "] Can't find general.in; exiting early.");
741  }
742 
743  GetPot infile("general.in");
744 
745  // Read in parameters from the input file
746  FEMParameters param(init.comm());
747  param.read(infile);
748 
749  // Skip higher-dimensional examples on a lower-dimensional libMesh build
750  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
751 
752  // Create a mesh with the given dimension, distributed
753  // across the default MPI communicator.
754  Mesh mesh(init.comm(), param.dimension);
755 
756  // And an object to refine it
757  UniquePtr<MeshRefinement> mesh_refinement =
758  build_mesh_refinement(mesh, param);
759 
760  // And an EquationSystems to run on it
761  EquationSystems equation_systems (mesh);
762 
763  libMesh::out << "Building mesh" << std::endl;
764 
765  if (!param.initial_timestep && param.run_simulation)
766  build_domain(mesh, param);
767 
768  libMesh::out << "Building system" << std::endl;
769 
770  FEMSystem & system = build_system(equation_systems, infile, param);
771 
772  set_system_parameters(system, param);
773 
774  libMesh::out << "Initializing systems" << std::endl;
775 
776  // Create a QoI object, we will need this to compute the QoI
777  {
778  CoupledSystemQoI qoi;
779  // Our QoI is computed on the side
780  qoi.assemble_qoi_sides = true;
781  system.attach_qoi(&qoi);
782  }
783 
784  equation_systems.init ();
785 
786  // Print information about the mesh and system to the screen.
787  mesh.print_info();
788  equation_systems.print_info();
789 
790  libMesh::out << "Starting adaptive loop" << std::endl << std::endl;
791 
792  // Count the number of steps used
793  unsigned int a_step = 0;
794 
795  // Get the linear solver object to set the preconditioner reuse flag
796  LinearSolver<Number> * linear_solver = system.get_linear_solver();
797 
798  // Adaptively solve the timestep
799  for (; a_step <= param.max_adaptivesteps; ++a_step)
800  {
801  // We have to ensure that we are refining to either an error
802  // tolerance or to a target number of elements, and, not to a
803  // non-zero tolerance and a target number of elements
804  // simultaneously
805  if (param.global_tolerance > 0 && param.nelem_target > 0)
806  {
807  libMesh::out<<"We cant refine to both a non-zero tolerance and a target number of elements, EXITING adaptive loop. "<<std::endl<<std::endl;
808  break;
809  }
810 
811  // We dont want to reuse the preconditioner for the forward solve
812  linear_solver->reuse_preconditioner(false);
813 
814  libMesh::out << "Adaptive step " << a_step << std::endl;
815 
816  libMesh::out << "Solving the forward problem" << std::endl;
817 
818  libMesh::out << "We have "
819  << mesh.n_active_elem()
820  << " active elements and "
821  << equation_systems.n_active_dofs()
822  << " active dofs."
823  << std::endl
824  << std::endl;
825 
826  // Solve the forward system
827  system.solve();
828 
829  // Write the output
830  write_output(equation_systems, 0, a_step, "primal", param);
831 
832  // Compute the QoI
833  system.assemble_qoi();
834 
835  // We just call a postprocess here to set the variable computed_QoI to the value computed by assemble_qoi
836  system.postprocess();
837 
838  // Get the value of the computed_QoI variable of the CoupledSystem class
839  Number QoI_0_computed = (dynamic_cast<CoupledSystem &>(system)).get_QoI_value();
840 
841  libMesh::out << "The boundary QoI is "
842  << std::setprecision(17)
843  << QoI_0_computed
844  << std::endl
845  << std::endl;
846 
847  // Hard coded assert to ensure that the actual numbers we are getting are what they should be
848  if (a_step == param.max_adaptivesteps)
849  libmesh_assert_less(std::abs(QoI_0_computed - 0.0833), 2.5e-4);
850 
851  // You dont need to compute error estimates and refine at the last
852  // adaptive step, only before that
853  if (a_step!=param.max_adaptivesteps)
854  {
855  NumericVector<Number> & primal_solution = (*system.solution);
856 
857  // Make sure we get the contributions to the adjoint RHS from the sides
858  system.assemble_qoi_sides = true;
859 
860  // We are about to solve the adjoint system, but before we
861  // do this we see the same preconditioner flag to reuse the
862  // preconditioner from the forward solver
863  linear_solver->reuse_preconditioner(param.reuse_preconditioner);
864 
865  // Solve the adjoint system
866  libMesh::out<< "Solving the adjoint problem" <<std::endl;
867  system.adjoint_solve();
868 
869  // Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
870  system.set_adjoint_already_solved(true);
871 
872  // To plot the adjoint solution, we swap it with the primal solution
873  // and use the write_output function
874  NumericVector<Number> & dual_solution = system.get_adjoint_solution();
875  primal_solution.swap(dual_solution);
876 
877  write_output(equation_systems, 0, a_step, "adjoint", param);
878 
879  // Swap back
880  primal_solution.swap(dual_solution);
881 
882  // We need the values of the parameters Pe from the
883  // system for the adjoint error estimate
884  Real Pe = (dynamic_cast<CoupledSystem &>(system)).get_Pe();
885 
886  // The total error is the sum: error = error_non_pressure +
887  // error_with_pressure + ...
888  // error_estimator_convection_diffusion_x +
889  // error_estimator_convection_diffusion_y
890  ErrorVector error;
891 
892  // We first construct the non-pressure contributions
893  ErrorVector error_non_pressure;
894 
895  // First we build the norm_type_vector_non_pressure vectors and
896  // weights_matrix_non_pressure matrix for the non-pressure term
897  // error contributions
898  std::vector<FEMNormType> primal_norm_type_vector_non_pressure;
899  primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
900  primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
901  primal_norm_type_vector_non_pressure.push_back(L2);
902  primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
903 
904  std::vector<FEMNormType> dual_norm_type_vector_non_pressure;
905  dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
906  dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
907  dual_norm_type_vector_non_pressure.push_back(L2);
908  dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
909 
910  std::vector<std::vector<Real>>
911  weights_matrix_non_pressure(system.n_vars(),
912  std::vector<Real>(system.n_vars(), 0.0));
913  weights_matrix_non_pressure[0][0] = 1.;
914  weights_matrix_non_pressure[1][1] = 1.;
915  weights_matrix_non_pressure[3][3] = 1./Pe;
916 
917  // We build the error estimator to estimate the contributions
918  // to the QoI error from the non pressure term
919  UniquePtr<ErrorEstimator> error_estimator_non_pressure =
920  build_error_estimator_component_wise (param,
921  weights_matrix_non_pressure,
922  primal_norm_type_vector_non_pressure,
923  dual_norm_type_vector_non_pressure);
924 
925  // Estimate the contributions to the QoI error from the non
926  // pressure terms
927  error_estimator_non_pressure->estimate_error(system, error_non_pressure);
928 
929  // Plot the estimated error from the non_pressure terms
930  write_error(equation_systems, error_non_pressure, 0, a_step, param, "_non_pressure");
931 
932  // Now for the pressure contributions
933  ErrorVector error_with_pressure;
934 
935  // Next we build the norm_type_vector_with_pressure vectors and
936  // weights_matrix_with_pressure matrix for the pressure term
937  // error contributions
938  std::vector<FEMNormType> primal_norm_type_vector_with_pressure;
939  primal_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
940  primal_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
941  primal_norm_type_vector_with_pressure.push_back(L2);
942  primal_norm_type_vector_with_pressure.push_back(L2);
943 
944  std::vector<FEMNormType> dual_norm_type_vector_with_pressure;
945  dual_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
946  dual_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
947  dual_norm_type_vector_with_pressure.push_back(L2);
948  dual_norm_type_vector_with_pressure.push_back(L2);
949 
950  std::vector<std::vector<Real>>
951  weights_matrix_with_pressure (system.n_vars(),
952  std::vector<Real>(system.n_vars(), 0.0));
953  weights_matrix_with_pressure[0][2] = 1.;
954  weights_matrix_with_pressure[1][2] = 1.;
955  weights_matrix_with_pressure[2][0] = 1.;
956  weights_matrix_with_pressure[2][1] = 1.;
957 
958  // We build the error estimator to estimate the contributions
959  // to the QoI error from the pressure term
960  UniquePtr<ErrorEstimator> error_estimator_with_pressure =
961  build_error_estimator_component_wise (param,
962  weights_matrix_with_pressure,
963  primal_norm_type_vector_with_pressure,
964  dual_norm_type_vector_with_pressure);
965 
966  // Estimate the contributions to the QoI error from the pressure terms
967  error_estimator_with_pressure->estimate_error(system, error_with_pressure);
968 
969  // Plot the error due to the pressure terms
970  write_error(equation_systems, error_with_pressure, 0, a_step, param, "_with_pressure");
971 
972  // Now for the convection diffusion term errors (in the x and y directions)
973 
974  ErrorVector error_convection_diffusion_x;
975 
976  // The norm type vectors and weights matrix for the convection_diffusion_x errors
977  std::vector<FEMNormType> primal_norm_type_vector_convection_diffusion_x;
978  primal_norm_type_vector_convection_diffusion_x.push_back(L2);
979  primal_norm_type_vector_convection_diffusion_x.push_back(L2);
980  primal_norm_type_vector_convection_diffusion_x.push_back(L2);
981  primal_norm_type_vector_convection_diffusion_x.push_back(H1_X_SEMINORM);
982 
983  std::vector<FEMNormType> dual_norm_type_vector_convection_diffusion_x;
984  dual_norm_type_vector_convection_diffusion_x.push_back(L2);
985  dual_norm_type_vector_convection_diffusion_x.push_back(L2);
986  dual_norm_type_vector_convection_diffusion_x.push_back(L2);
987  // Note that we need the error of the dual concentration in L2
988  dual_norm_type_vector_convection_diffusion_x.push_back(L2);
989 
990  std::vector<std::vector<Real>>
991  weights_matrix_convection_diffusion_x (system.n_vars(),
992  std::vector<Real>(system.n_vars(), 0.0));
993  weights_matrix_convection_diffusion_x[0][3] = 1.;
994  weights_matrix_convection_diffusion_x[3][3] = 1.;
995 
996  // We will also have to build and pass the weight functions to the weighted patch recovery estimators
997 
998  // We pass the identity function as weights to error entries that the above matrix will scale to 0.
999  ConstFEMFunction<Number> identity(1);
1000 
1001  // Declare object of class CoupledFEMFunctionsx, the definition of the function contains the weight
1002  // to be applied to the relevant terms
1003  // For ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2} term, returns C,1_h
1004  CoupledFEMFunctionsx convdiffx0(system, 0);
1005  // For ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} term, returns (u_1)_h
1006  CoupledFEMFunctionsx convdiffx3(system, 3);
1007 
1008  // Make a vector of pointers to these objects
1009  std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_x;
1010  coupled_system_weight_functions_x.push_back(&convdiffx0);
1011  coupled_system_weight_functions_x.push_back(&identity);
1012  coupled_system_weight_functions_x.push_back(&identity);
1013  coupled_system_weight_functions_x.push_back(&convdiffx3);
1014 
1015  // Build the error estimator to estimate the contributions
1016  // to the QoI error from the convection diffusion x term
1017  UniquePtr<ErrorEstimator> error_estimator_convection_diffusion_x =
1018  build_weighted_error_estimator_component_wise (param,
1019  weights_matrix_convection_diffusion_x,
1020  primal_norm_type_vector_convection_diffusion_x,
1021  dual_norm_type_vector_convection_diffusion_x,
1022  coupled_system_weight_functions_x);
1023 
1024  // Estimate the contributions to the QoI error from the
1025  // convection diffusion x term
1026  error_estimator_convection_diffusion_x->estimate_error(system, error_convection_diffusion_x);
1027 
1028  // Plot this error
1029  write_error (equation_systems,
1030  error_convection_diffusion_x,
1031  0,
1032  a_step,
1033  param,
1034  "_convection_diffusion_x");
1035 
1036  // Now for the y direction terms
1037  ErrorVector error_convection_diffusion_y;
1038 
1039  // The norm type vectors and weights matrix for the convection_diffusion_x errors
1040  std::vector<FEMNormType> primal_norm_type_vector_convection_diffusion_y;
1041  primal_norm_type_vector_convection_diffusion_y.push_back(L2);
1042  primal_norm_type_vector_convection_diffusion_y.push_back(L2);
1043  primal_norm_type_vector_convection_diffusion_y.push_back(L2);
1044  primal_norm_type_vector_convection_diffusion_y.push_back(H1_Y_SEMINORM);
1045 
1046  std::vector<FEMNormType> dual_norm_type_vector_convection_diffusion_y;
1047  dual_norm_type_vector_convection_diffusion_y.push_back(L2);
1048  dual_norm_type_vector_convection_diffusion_y.push_back(L2);
1049  dual_norm_type_vector_convection_diffusion_y.push_back(L2);
1050  // Note that we need the error of the dual concentration in L2
1051  dual_norm_type_vector_convection_diffusion_y.push_back(L2);
1052 
1053  std::vector<std::vector<Real>>
1054  weights_matrix_convection_diffusion_y (system.n_vars(),
1055  std::vector<Real>(system.n_vars(), 0.0));
1056  weights_matrix_convection_diffusion_y[1][3] = 1.;
1057  weights_matrix_convection_diffusion_y[3][3] = 1.;
1058 
1059  // For ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2} term, returns C,2_h
1060  CoupledFEMFunctionsy convdiffy1(system, 1);
1061  // For ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} term, returns (u_2)_h
1062  CoupledFEMFunctionsy convdiffy3(system, 3);
1063 
1064  // Make a vector of pointers to these objects
1065  std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_y;
1066  coupled_system_weight_functions_y.push_back(&identity);
1067  coupled_system_weight_functions_y.push_back(&convdiffy1);
1068  coupled_system_weight_functions_y.push_back(&identity);
1069  coupled_system_weight_functions_y.push_back(&convdiffy3);
1070 
1071  // Build the error estimator to estimate the contributions
1072  // to the QoI error from the convection diffusion y term
1073  UniquePtr<ErrorEstimator> error_estimator_convection_diffusion_y =
1074  build_weighted_error_estimator_component_wise (param,
1075  weights_matrix_convection_diffusion_y,
1076  primal_norm_type_vector_convection_diffusion_y,
1077  dual_norm_type_vector_convection_diffusion_y,
1078  coupled_system_weight_functions_y);
1079 
1080  // Estimate the contributions to the QoI error from the
1081  // convection diffusion y terms
1082  error_estimator_convection_diffusion_y->estimate_error(system, error_convection_diffusion_y);
1083 
1084  // Plot this error
1085  write_error(equation_systems, error_convection_diffusion_y, 0, a_step, param, "_convection_diffusion_y");
1086 
1087  if (param.indicator_type == "adjoint_residual")
1088  {
1089  error.resize(error_non_pressure.size());
1090 
1091  // Now combine the contribs from the pressure and non
1092  // pressure terms to get the complete estimate
1093  for (std::size_t i = 0; i < error.size(); i++)
1094  error[i] = error_non_pressure[i] + error_with_pressure[i] + error_convection_diffusion_x[i] + error_convection_diffusion_y[i];
1095  }
1096  else
1097  {
1098  libMesh::out << "Using Kelly Estimator" << std::endl;
1099 
1100  // Build the Kelly error estimator
1101  UniquePtr<ErrorEstimator> error_estimator = build_error_estimator(param);
1102 
1103  // Estimate the error
1104  error_estimator->estimate_error(system, error);
1105  }
1106 
1107  write_error(equation_systems, error, 0, a_step, param, "_total");
1108 
1109  // We have to refine either based on reaching an error
1110  // tolerance or a number of elements target, which should be
1111  // verified above
1112 
1113  // Otherwise we flag elements by error tolerance or nelem
1114  // target
1115 
1116  if (param.refine_uniformly)
1117  {
1118  mesh_refinement->uniformly_refine(1);
1119  }
1120  else if (param.global_tolerance >= 0. && param.nelem_target == 0.) // Refine based on reaching an error tolerance
1121  {
1122  mesh_refinement->flag_elements_by_error_tolerance (error);
1123 
1124  // Carry out the adaptive mesh refinement/coarsening
1125  mesh_refinement->refine_and_coarsen_elements();
1126  }
1127  else // Refine based on reaching a target number of elements
1128  {
1129  //If we have reached the desired number of elements, we
1130  //dont do any more adaptive steps
1131  if (mesh.n_active_elem() >= param.nelem_target)
1132  {
1133  libMesh::out << "We reached the target number of elements." << std::endl <<std::endl;
1134  break;
1135  }
1136 
1137  mesh_refinement->flag_elements_by_nelem_target (error);
1138 
1139  // Carry out the adaptive mesh refinement/coarsening
1140  mesh_refinement->refine_and_coarsen_elements();
1141  }
1142 
1143  equation_systems.reinit();
1144  }
1145  }
1146 
1147  write_output_footers(param);
1148 
1149 #endif // #ifndef LIBMESH_ENABLE_AMR
1150 
1151  // All done.
1152  return 0;
1153 
1154 }
bool continue_after_max_iterations
Defaults to true, telling the DiffSolver to continue rather than exit when a solve has reached its ma...
Definition: diff_solver.h:174
unsigned int nelem_target
Definition: femparameters.h:56
libMesh::Real timesolver_upper_tolerance
Definition: femparameters.h:38
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1545
bool print_solution_norms
void start_output(unsigned int timesteps, std::string filename, std::string varname)
Definition: output.C:11
This class implements a goal oriented error indicator, by weighting residual-based estimates from the...
libMesh::Real deltat
Definition: femparameters.h:38
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
Real minimum_linear_tolerance
The tolerance for linear solves is kept above this minimum.
Definition: diff_solver.h:215
This is the EquationSystems class.
unsigned int initial_timestep
Definition: femparameters.h:34
bool print_element_residuals
libMesh::Real initial_linear_tolerance
unsigned int dimension
Definition: femparameters.h:45
void write_output_headers(FEMParameters &param)
Definition: adjoints_ex3.C:250
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:203
bool reuse_preconditioner
void write(const std::string &name, const XdrMODE, const unsigned int write_flags=(WRITE_DATA), bool partition_agnostic=true) const
Write the systems to disk using the XDR data format.
libMesh::Real timesolver_tolerance
Definition: femparameters.h:38
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time)
Writes out the solution at a specific timestep.
Definition: exodusII_io.C:773
void set_adjoint_already_solved(bool setting)
Setter for the adjoint_already_solved boolean.
Definition: system.h:378
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
SystemNorm component_norm
Error calculations are done in this norm, DISCRETE_L2 by default.
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:162
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Definition: diff_system.h:329
Real & absolute_global_tolerance()
If absolute_global_tolerance is set to a nonzero value, methods like flag_elements_by_global_toleranc...
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
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
int extra_quadrature_order
A member int that can be employed to indicate increased or reduced quadrature order.
Definition: system.h:1508
Real theta
The value for the theta method to employ: 1.0 corresponds to backwards Euler, 0.0 corresponds to forw...
Definition: euler_solver.h:100
std::string timesolver_core
Definition: femparameters.h:37
unsigned int max_nonlinear_iterations
The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_...
Definition: diff_solver.h:156
Real absolute_residual_tolerance
The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolera...
Definition: diff_solver.h:191
bool require_residual_reduction
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
MeshBase & mesh
bool print_residual_norms
libMesh::Real timesolver_maxgrowth
Definition: femparameters.h:38
Real target_tolerance
This tolerance is the target relative error between an exact time integration and a single time step ...
std::string numbered_filename(unsigned int t_step, unsigned int a_step, std::string solution_type, std::string type, std::string extension, FEMParameters &param)
Definition: adjoints_ex3.C:119
Real & refine_fraction()
The refine_fraction sets either a desired target or a desired maximum number of elements to flag for ...
libMesh::Real relative_step_tolerance
std::string indicator_type
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
Definition: system_norm.h:43
UniquePtr< ErrorEstimator > & primal_error_estimator()
Access to the "subestimator" (default: PatchRecovery) to use on the primal/forward solution...
void read(const std::string &name, const XdrMODE, const unsigned int read_flags=(READ_HEADER|READ_DATA), bool partition_agnostic=true)
Read & initialize the systems from disk using the XDR data format.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
Real initial_linear_tolerance
Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the ...
Definition: diff_solver.h:210
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
libMesh::Real refine_fraction
Definition: femparameters.h:58
void set_system_parameters(LaplaceSystem &system, FEMParameters &param)
Definition: adjoints_ex1.C:161
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
Definition: diff_system.C:175
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:334
This class provides a specific system class.
Definition: fem_system.h:53
bool & coarsen_by_parents()
If coarsen_by_parents is true, complete groups of sibling elements (elements with the same parent) wi...
This class defines a solver which uses the default libMesh linear solver in a quasiNewton method to h...
Definition: newton_solver.h:46
This is the MeshBase class.
Definition: mesh_base.h:68
void set_off_diagonal_weight(unsigned int i, unsigned int j, Real w)
Sets the weight corresponding to the norm from the variable pair v1(var1) coming from v2(var2)...
Definition: system_norm.h:316
bool print_element_residuals
Set print_element_residuals to true to print each R_elem contribution.
Definition: diff_system.h:344
libmesh_assert(j)
dof_id_type & nelem_target()
If nelem_target is set to a nonzero value, methods like flag_elements_by_nelem_target() will attempt ...
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscDiffSolver & solver
bool print_element_solutions
Set print_element_solutions to true to print each U_elem input.
Definition: diff_system.h:339
SolverPackage default_solver_package()
Definition: libmesh.C:995
This class defines a solver which uses a PETSc SNES context to handle a DifferentiableSystem.
libMesh::Real coarsen_threshold
Definition: femparameters.h:58
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
bool write_gmv_error
Definition: femparameters.h:65
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
This class defines a theta-method (defaulting to Backward Euler with theta = 1.0) solver to handle ti...
Definition: euler2_solver.h:48
This is the MeshRefinement class.
unsigned int max_linear_iterations
bool print_residual_norms
Set print_residual_norms to true to print |F| whenever it is assembled.
Definition: diff_system.h:319
void set_weight(unsigned int var, Real w)
Sets the weight corresponding to the norm in variable var.
Definition: system_norm.h:301
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:148
Real upper_tolerance
This tolerance is the maximum relative error between an exact time integration and a single time step...
libMesh::Real minimum_linear_tolerance
virtual void reinit()
Reinitialize all the systems.
libMesh::Real verify_analytic_jacobians
Real & coarsen_threshold()
The coarsen_threshold provides hysteresis in AMR/C strategies.
bool print_residuals
Set print_residuals to true to print F whenever it is assembled.
Definition: diff_system.h:324
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
unsigned int max_adaptivesteps
Definition: femparameters.h:59
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void plot_error(const std::string &filename, const MeshBase &mesh) const
Plots a data file, of a type determined by looking at the file extension in filename, of the error values on the active elements of mesh.
Definition: error_vector.C:208
int extra_quadrature_order
libMesh::Real global_tolerance
Definition: femparameters.h:57
virtual void write(const std::string &name)=0
Real theta
The value for the theta method to employ: 1.0 corresponds to backwards Euler, 0.0 corresponds to forw...
void attach_qoi(DifferentiableQoI *qoi_in)
Attach external QoI object.
Definition: diff_system.h:212
Real linear_tolerance_multiplier
The tolerance for linear solves is kept below this multiplier (which defaults to 1e-3) times the norm...
libMesh::Real absolute_residual_tolerance
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
void write_output_footers(FEMParameters &param)
Definition: adjoints_ex3.C:309
unsigned int n_systems() const
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
This class defines a theta-method Euler (defaulting to Backward Euler with theta = 1...
Definition: euler_solver.h:43
This class implements the Kelly error indicator which is based on the flux jumps between elements...
bool verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
Definition: diff_solver.h:168
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
bool print_element_solutions
void build_domain(MeshBase &mesh, FEMParameters &param)
Definition: domain.C:10
This is a generic class that defines a solver to handle time integration of DifferentiableSystems.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This class wraps another UnsteadySolver derived class, and compares the results of timestepping with ...
bool continue_after_backtrack_failure
Defaults to false, telling the DiffSolver to throw an error when the backtracking scheme fails to fin...
Definition: diff_solver.h:180
void write_output_solvedata(EquationSystems &es, unsigned int a_step, unsigned int newton_steps, unsigned int krylov_steps, unsigned int tv_sec, unsigned int tv_usec)
Definition: adjoints_ex3.C:273
virtual void reuse_preconditioner(bool)
Set the same_preconditioner flag, which indicates if we reuse the same preconditioner for subsequent ...
bool print_jacobian_norms
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
Real minsteplength
If the quasi-Newton step length must be reduced to below this factor to give a residual reduction...
bool assemble_qoi_sides
If assemble_qoi_sides is true (it is false by default), the assembly loop for a quantity of interest ...
Definition: diff_qoi.h:82
This class implements the Patch Recovery error indicator.
FEMSystem & build_system(EquationSystems &es, GetPot &, FEMParameters &)
Definition: mysystems.h:7
OStreamProxy out
processor_id_type global_processor_id()
Definition: libmesh_base.h:114
void write_error(EquationSystems &es, ErrorVector &error, unsigned int t_number, unsigned int a_number, FEMParameters &param, std::string error_type) void write_error(EquationSystems &
libMesh::Real min_step_length
libMesh::Real timesolver_theta
Definition: femparameters.h:38
std::vector< FEMFunctionBase< Number > * > weight_functions
Vector of fem function base pointers, the user will fill this in with pointers to the appropriate wei...
void read(GetPot &input, const std::vector< std::string > *other_variable_names=libmesh_nullptr)
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:989
UniquePtr< ErrorEstimator > build_error_estimator(FEMParameters &param, QoISet &qois)
Definition: adjoints_ex1.C:239
const MeshBase & get_mesh() const
unsigned int max_nonlinear_iterations
This class implements writing meshes in the Tecplot format.
Definition: tecplot_io.h:43
unsigned int rank() const
Definition: parallel.h:724
UniquePtr< MeshRefinement > build_mesh_refinement(MeshBase &mesh, FEMParameters &param)
Definition: adjoints_ex1.C:216
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
bool require_residual_reduction
If this is set to true, the solver is forced to test the residual after each Newton step...
Definition: newton_solver.h:98
FEMFunction that returns a single value, regardless of the time and location inputs.
void write_output(EquationSystems &es, unsigned int t_step, unsigned int a_step, std::string solution_type, FEMParameters &param)
Definition: adjoints_ex3.C:151
bool time_solver_quiet
virtual void init()
Initialize all the systems.
std::size_t n_active_dofs() const
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
UniquePtr< ErrorEstimator > & dual_error_estimator()
Access to the "subestimator" (default: PatchRecovery) to use on the dual/adjoint solution.
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
UniquePtr< UnsteadySolver > core_time_solver
This object is used to take timesteps.
std::vector< libMesh::FEMNormType > timesolver_norm
Definition: femparameters.h:41
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
virtual void assemble_qoi(const QoISet &indices=QoISet()) libmesh_override
Runs a qoi assembly loop over all elements, and if assemble_qoi_sides is true over all sides...
Definition: fem_system.C:1127
libMesh::Real linear_tolerance_multiplier
This class implements the Patch Recovery error indicator.
virtual void renumber_nodes_and_elements()=0
After partitioning a mesh it is useful to renumber the nodes and elements so that they lie in contigu...
bool write_tecplot_error
Definition: femparameters.h:65
Real relative_residual_tolerance
Definition: diff_solver.h:192
libMesh::Real coarsen_fraction
Definition: femparameters.h:58
virtual void postprocess() libmesh_override
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: fem_system.C:1107
processor_id_type processor_id() const
int main(int argc, char **argv)
libMesh::Real steadystate_tolerance
Definition: femparameters.h:38
Real max_growth
Do not allow the adaptive time solver to select a new deltat greater than max_growth times the old de...
unsigned int deltat_reductions
Definition: femparameters.h:36