libMesh
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
20 // Adjoints Example 5 - SolutionHistory, General Localized Vectors and Unsteady Adjoints.
21 // \author Vikram Garg
22 // \date 2013
23 //
24 // This example showcases three new capabilities in libMesh. The
25 // primary motivation is adjoint sensitivity analysis for unsteady
26 // problems. The PDE we are interested in is the simple 2-d heat
27 // equation:
28 // partial(T)/partial(t) - K Laplacian(T) = 0
29 // with initial condition:
30 // T(x,y;0) = sin(pi*x) sin(pi*y) and boundary conditions:
31 // T(boundary;t) = 0
32
33 // For these initial and boundary conditions, the exact solution
34 // u = exp(-K pi^2 t) * sin(pi*x) * sin(pi*y)
35
36 // We specify our Quantity of Interest (QoI) as
37 // Q(u) = int_{domain} u(x,y;1) sin(pi*x) sin(pi*y) dx dy, and
38 // are interested in computing the sensitivity dQ/dK
39
40 // For this QoI, the continuous adjoint problem reads,
41 // -partial(z)/partial(t) - K Laplacian(z) = 0
42 // with initial condition:
43 // T(x,y;1) = sin(pi*x) sin(pi*y)
44 // and boundary condition:
45 // T(boundary;t) = 0
46
47 // which has the exact solution,
48 // z = exp(-K pi^2 (1 - t)) * sin(pi*x) * sin(pi*y)
49 // which is the mirror image in time of the forward solution
50
51 // For an adjoint consistent space-time formulation, the discrete
52 // adjoint can be obtained by marching backwards from the adjoint
53 // initial condition and solving the transpose of the discrete primal
54 // problem at the last nonlinear solve of the corresponding primal
55 // timestep. This necessitates the storage of the primal solution at
56 // all timesteps, which is accomplished here using a
57 // MemorySolutionHistory object. As the name suggests, this object
58 // simply stores the primal solution (and other vectors we may choose
59 // to save), so that we can retrieve them later, whenever necessary.
60
61 // The discrete adjoint system for implicit time steppers requires the
62 // localization of vectors other than system.solution, which is
63 // accomplished using the localize_vectors method. In this particular
64 // example, we use the localized adjoint solution to assemble the
65 // residual contribution for the current adjoint timestep from the last
66 // computed adjoint timestep.
67
68 // Finally, The adjoint_advance_timestep method, the backwards time
69 // analog of advance_timestep prepares the time solver for solving the
70 // adjoint system, while the retrieve_timestep method retrieves the
71 // saved solutions at the current system.time, so that the adjoint
72 // sensitivity contribution for the current time can be computed.
73
74 // Local includes
75 #include "initial.h"
77 #include "femparameters.h"
78 #include "heatsystem.h"
79
80 // Libmesh includes
81 #include "libmesh/equation_systems.h"
82 #include "libmesh/dirichlet_boundaries.h"
83 #include "libmesh/system_norm.h"
84 #include "libmesh/numeric_vector.h"
85
86 #include "libmesh/mesh.h"
87 #include "libmesh/mesh_generation.h"
88 #include "libmesh/mesh_refinement.h"
89
90 #include "libmesh/petsc_diff_solver.h"
92 #include "libmesh/euler_solver.h"
93 #include "libmesh/euler2_solver.h"
94 #include "libmesh/newton_solver.h"
95 #include "libmesh/twostep_time_solver.h"
96
97 #include "libmesh/getpot.h"
98 #include "libmesh/tecplot_io.h"
99 #include "libmesh/gmv_io.h"
100 #include "libmesh/exodusII_io.h"
101
102 // SolutionHistory Includes
103 #include "libmesh/solution_history.h"
104 #include "libmesh/memory_solution_history.h"
105
106 // C++ includes
107 #include <iostream>
108 #include <sys/time.h>
109 #include <iomanip>
110
111 void write_output(EquationSystems & es,
112  unsigned int t_step, // The current time step count
113  std::string solution_type, // primal or adjoint solve
114  FEMParameters & param)
115 {
116  // Ignore parameters when there are no output formats available.
117  libmesh_ignore(es);
118  libmesh_ignore(t_step);
119  libmesh_ignore(solution_type);
120  libmesh_ignore(param);
121
122 #ifdef LIBMESH_HAVE_GMV
123  if (param.output_gmv)
124  {
125  MeshBase & mesh = es.get_mesh();
126
127  std::ostringstream file_name_gmv;
128  file_name_gmv << solution_type
129  << ".out.gmv."
130  << std::setw(2)
131  << std::setfill('0')
132  << std::right
133  << t_step;
134
135  GMVIO(mesh).write_equation_systems(file_name_gmv.str(), es);
136  }
137 #endif
138
139 #ifdef LIBMESH_HAVE_EXODUS_API
140  if (param.output_exodus)
141  {
142  MeshBase & mesh = es.get_mesh();
143
144  // We write out one file per timestep. The files are named in
145  // the following way:
146  // foo.e
147  // foo.e-s002
148  // foo.e-s003
149  // ...
150  // so that, if you open the first one with Paraview, it actually
151  // opens the entire sequence of adapted files.
152  std::ostringstream file_name_exodus;
153
154  file_name_exodus << solution_type << ".e";
155  if (t_step > 0)
156  file_name_exodus << "-s"
157  << std::setw(3)
158  << std::setfill('0')
159  << std::right
160  << t_step + 1;
161
162  // TODO: Get the current time from the System...
163  ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
164  es,
165  1,
166  /*time=*/t_step + 1);
167  }
168 #endif
169 }
170
172  FEMParameters & param)
173 {
174  // Use the prescribed FE type
175  system.fe_family() = param.fe_family[0];
176  system.fe_order() = param.fe_order[0];
177
178  // Use analytical jacobians?
179  system.analytic_jacobians() = param.analytic_jacobians;
180
181  // Verify analytic jacobians against numerical ones?
184
185  // More desperate debugging options
187  system.print_solutions = param.print_solutions;
189  system.print_residuals = param.print_residuals;
191  system.print_jacobians = param.print_jacobians;
194
195  // Solve this as a time-dependent or steady system
196  if (param.transient)
197  {
199  if (param.timesolver_core == "euler")
200  {
201  EulerSolver *eulersolver =
202  new EulerSolver(system);
203
204  eulersolver->theta = param.timesolver_theta;
205  innersolver = eulersolver;
206  }
207  else
208  libmesh_error_msg("This example (and unsteady adjoints in libMesh) only support Backward Euler and explicit methods.");
209
210  system.time_solver =
211  UniquePtr<TimeSolver>(innersolver);
212  }
213  else
214  system.time_solver =
216
217  // The Memory Solution History object we will set the system SolutionHistory object to
218  MemorySolutionHistory heatsystem_solution_history(system);
219  system.time_solver->set_solution_history(heatsystem_solution_history);
220
221  system.time_solver->reduce_deltat_on_diffsolver_failure =
222  param.deltat_reductions;
223  system.time_solver->quiet = param.time_solver_quiet;
224
225  // Create any Dirichlet boundary conditions
226  typedef
227  std::map<boundary_id_type, FunctionBase<Number> *>::
228  const_iterator Iter;
229
230  for (Iter i = param.dirichlet_conditions.begin();
231  i != param.dirichlet_conditions.end(); ++i)
232  {
233  boundary_id_type b = i->first;
234  FunctionBase<Number> *f = i->second;
235  std::set<boundary_id_type> bdys; bdys.insert(b);
236
239  f));
240
241  libMesh::out << "Added Dirichlet boundary " << b << " for variables ";
242  for (std::size_t vi=0; vi != param.dirichlet_condition_variables[b].size(); ++vi)
244  libMesh::out << std::endl;
245  }
246
247  // Set the time stepping options
248  system.deltat = param.deltat;
249
250  // And the integration options
252
253  // And the nonlinear solver options
254  if (param.use_petsc_snes)
255  {
256 #ifdef LIBMESH_HAVE_PETSC
257  PetscDiffSolver *solver = new PetscDiffSolver(system);
258  system.time_solver->diff_solver() = UniquePtr<DiffSolver>(solver);
259 #else
260  libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
261 #endif
262  }
263  else
264  {
265  NewtonSolver *solver = new NewtonSolver(system);
266  system.time_solver->diff_solver() = UniquePtr<DiffSolver>(solver);
267
268  solver->quiet = param.solver_quiet;
269  solver->verbose = param.solver_verbose;
270  solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
271  solver->minsteplength = param.min_step_length;
272  solver->relative_step_tolerance = param.relative_step_tolerance;
273  solver->relative_residual_tolerance = param.relative_residual_tolerance;
274  solver->require_residual_reduction = param.require_residual_reduction;
275  solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
276  if (system.time_solver->reduce_deltat_on_diffsolver_failure)
277  {
278  solver->continue_after_max_iterations = true;
279  solver->continue_after_backtrack_failure = true;
280  }
281
282  // And the linear solver options
283  solver->max_linear_iterations = param.max_linear_iterations;
284  solver->initial_linear_tolerance = param.initial_linear_tolerance;
285  solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
286  }
287 }
288
289 // The main program.
290 int main (int argc, char ** argv)
291 {
292  // Skip adaptive examples on a non-adaptive libMesh build
293 #ifndef LIBMESH_ENABLE_AMR
294  libmesh_ignore(argc);
295  libmesh_ignore(argv);
296  libmesh_example_requires(false, "--enable-amr");
297 #else
298  // Skip this 2D example if libMesh was compiled as 1D-only.
299  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
300
301  // Initialize libMesh.
302  LibMeshInit init (argc, argv);
303
304  // This doesn't converge with Trilinos for some reason...
305  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
306
307  libMesh::out << "Started " << argv[0] << std::endl;
308
309  // Make sure the general input file exists, and parse it
310  {
311  std::ifstream i("general.in");
312  if (!i)
313  libmesh_error_msg('[' << init.comm().rank() << "] Can't find general.in; exiting early.");
314  }
315  GetPot infile("general.in");
316
317  // Read in parameters from the input file
318  FEMParameters param(init.comm());
320
321  // Create a mesh with the given dimension, distributed
322  // across the default MPI communicator.
323  Mesh mesh(init.comm(), param.dimension);
324
325  // And an object to refine it
326  UniquePtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
327
328  // And an EquationSystems to run on it
329  EquationSystems equation_systems (mesh);
330
331  libMesh::out << "Building mesh" << std::endl;
332
333  // Build a unit square
334  ElemType elemtype;
335
336  if (param.elementtype == "tri" ||
337  param.elementtype == "unstructured")
338  elemtype = TRI3;
339  else
340  elemtype = QUAD4;
341
342  MeshTools::Generation::build_square (mesh, param.coarsegridx, param.coarsegridy,
343  param.domain_xmin, param.domain_xmin + param.domain_edge_width,
344  param.domain_ymin, param.domain_ymin + param.domain_edge_length,
345  elemtype);
346
347  libMesh::out << "Building system" << std::endl;
348
349  HeatSystem & system = equation_systems.add_system<HeatSystem> ("HeatSystem");
350
351  set_system_parameters(system, param);
352
353  libMesh::out << "Initializing systems" << std::endl;
354
355  // Initialize the system
356  equation_systems.init ();
357
358  // Refine the grid again if requested
359  for (unsigned int i=0; i != param.extrarefinements; ++i)
360  {
361  mesh_refinement->uniformly_refine(1);
362  equation_systems.reinit();
363  }
364
365  libMesh::out << "Setting primal initial conditions" << std::endl;
366
368
370  equation_systems.parameters);
371
372  // Output the H1 norm of the initial conditions
373  libMesh::out << "|U("
374  << system.time
375  << ")|= "
376  << system.calculate_norm(*system.solution, 0, H1)
377  << std::endl
378  << std::endl;
379
380  // Add an adjoint vector, this will be computed after the forward
381  // time stepping is complete
382  //
383  // Tell the library not to save adjoint solutions during the forward
384  // solve
385  //
386  // Tell the library not to project this vector, and hence, memory
387  // solution history to not save it.
388  //
389  // Make this vector ghosted so we can localize it to each element
390  // later.
391  const std::string & adjoint_solution_name = "adjoint_solution0";
393
394  // Close up any resources initial.C needed
396
397  // Plot the initial conditions
398  write_output(equation_systems, 0, "primal", param);
399
400  // Print information about the mesh and system to the screen.
401  mesh.print_info();
402  equation_systems.print_info();
403
404  // In optimized mode we catch any solver errors, so that we can
405  // write the proper footers before closing. In debug mode we just
406  // let the exception throw so that gdb can grab it.
407 #ifdef NDEBUG
408  try
409  {
410 #endif
411  // Now we begin the timestep loop to compute the time-accurate
412  // solution of the equations.
413  for (unsigned int t_step=param.initial_timestep;
414  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
415  {
416  // A pretty update message
417  libMesh::out << " Solving time step "
418  << t_step
419  << ", time = "
420  << system.time
421  << std::endl;
422
423  // Solve the forward problem at time t, to obtain the solution at time t + dt
424  system.solve();
425
426  // Output the H1 norm of the computed solution
427  libMesh::out << "|U("
428  << system.time + system.deltat
429  << ")|= "
430  << system.calculate_norm(*system.solution, 0, H1)
431  << std::endl;
432
433  // Advance to the next timestep in a transient problem
434  libMesh::out << "Advancing timestep" << std::endl << std::endl;
436
437  // Write out this timestep
438  write_output(equation_systems, t_step+1, "primal", param);
439  }
440  // End timestep loop
441
443
444  // Now we will solve the backwards in time adjoint problem
445  libMesh::out << std::endl << "Solving the adjoint problem" << std::endl;
446
447  // We need to tell the library that it needs to project the adjoint, so
448  // MemorySolutionHistory knows it has to save it
449
450  // Tell the library to project the adjoint vector, and hence, memory solution history to
451  // save it
453
454  libMesh::out << "Setting adjoint initial conditions Z("
455  << system.time
456  << ")"
457  <<std::endl;
458
459  // Need to call adjoint_advance_timestep once for the initial condition setup
460  libMesh::out<<"Retrieving solutions at time t="<<system.time<<std::endl;
462
463  // Output the H1 norm of the retrieved solutions (u^i and u^i+1)
464  libMesh::out << "|U("
465  << system.time + system.deltat
466  << ")|= "
467  << system.calculate_norm(*system.solution, 0, H1)
468  << std::endl;
469
470  libMesh::out << "|U("
471  << system.time
472  << ")|= "
473  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
474  << std::endl;
475
476  // The first thing we have to do is to apply the adjoint initial
477  // condition. The user should supply these. Here they are specified
481  equation_systems.parameters,
483
484  // Since we have specified an adjoint solution for the current
485  // time (T), set the adjoint_already_solved boolean to true, so
486  // we dont solve unnecessarily in the adjoint sensitivity method
488
489  libMesh::out << "|Z("
490  << system.time
491  << ")|= "
492  << system.calculate_norm(system.get_adjoint_solution(), 0, H1)
493  << std::endl
494  << std::endl;
495
496  write_output(equation_systems, param.n_timesteps, "dual", param);
497
498  // Now that the adjoint initial condition is set, we will start the
499  // backwards in time adjoint integration
500
501  // For loop stepping backwards in time
502  for (unsigned int t_step=param.initial_timestep;
503  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
504  {
505  //A pretty update message
506  libMesh::out << " Solving adjoint time step "
507  << t_step
508  << ", time = "
509  << system.time
510  << std::endl;
511
512  // The adjoint_advance_timestep function calls the retrieve
513  // function of the memory_solution_history class via the
514  // memory_solution_history object we declared earlier. The
515  // retrieve function sets the system primal vectors to their
516  // values at the current timestep.
517  libMesh::out << "Retrieving solutions at time t=" << system.time << std::endl;
519
520  // Output the H1 norm of the retrieved solution
521  libMesh::out << "|U("
522  << system.time + system.deltat
523  << ")|= "
524  << system.calculate_norm(*system.solution, 0, H1)
525  << std::endl;
526
527  libMesh::out << "|U("
528  << system.time
529  << ")|= "
530  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
531  << std::endl;
532
534
536
537  // Now that we have solved the adjoint, set the
538  // adjoint_already_solved boolean to true, so we dont solve
539  // unnecessarily in the error estimator
541
542  libMesh::out << "|Z("
543  << system.time
544  << ")|= "
545  << system.calculate_norm(system.get_adjoint_solution(), 0, H1)
546  << std::endl
547  << std::endl;
548
549  // Get a pointer to the primal solution vector
550  NumericVector<Number> & primal_solution = *system.solution;
551
552  // Get a pointer to the solution vector of the adjoint problem for QoI 0
553  NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
554
555  // Swap the primal and dual solutions so we can write out the adjoint solution
556  primal_solution.swap(dual_solution_0);
557
558  write_output(equation_systems, param.n_timesteps - (t_step + 1), "dual", param);
559
560  // Swap back
561  primal_solution.swap(dual_solution_0);
562  }
563  // End adjoint timestep loop
564
565  // Now that we have computed both the primal and adjoint solutions, we compute the sensitivities to the parameter p
566  // dQ/dp = partialQ/partialp - partialR/partialp
567  // partialQ/partialp = (Q(p+dp) - Q(p-dp))/(2*dp), this is not supported by the library yet
568  // partialR/partialp = (R(u,z;p+dp) - R(u,z;p-dp))/(2*dp), where
569  // R(u,z;p+dp) = int_{0}^{T} f(z;p+dp) - <partialu/partialt, z>(p+dp) - <g(u),z>(p+dp)
570  // To do this we need to step forward in time, and compute the perturbed R at each time step and accumulate it
571  // Then once all time steps are over, we can compute (R(u,z;p+dp) - R(u,z;p-dp))/(2*dp)
572
573  // Now we begin the timestep loop to compute the time-accurate
574  // adjoint sensitivities
575  for (unsigned int t_step=param.initial_timestep;
576  t_step != param.initial_timestep + param.n_timesteps; ++t_step)
577  {
578  // A pretty update message
579  libMesh::out << "Retrieving "
580  << t_step
581  << ", time = "
582  << system.time
583  << std::endl;
584
585  // Retrieve the primal and adjoint solutions at the current timestep
586  system.time_solver->retrieve_timestep();
587
588  libMesh::out << "|U("
589  << system.time + system.deltat
590  << ")|= "
591  << system.calculate_norm(*system.solution, 0, H1)
592  << std::endl;
593
594  libMesh::out << "|U("
595  << system.time
596  << ")|= "
597  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
598  << std::endl;
599
600  libMesh::out << "|Z("
601  << system.time
602  << ")|= "
603  << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
604  << std::endl
605  << std::endl;
606
607  // Call the postprocess function which we have overloaded to compute
608  // accumulate the perturbed residuals
609  dynamic_cast<HeatSystem &>(system).perturb_accumulate_residuals(dynamic_cast<HeatSystem &>(system).get_parameter_vector());
610
611  // Move the system time forward (retrieve_timestep does not do this)
612  system.time += system.deltat;
613  }
614
615  // A pretty update message
616  libMesh::out << "Retrieving final time = "
617  << system.time
618  << std::endl;
619
620  // Retrieve the primal and adjoint solutions at the current timestep
621  system.time_solver->retrieve_timestep();
622
623  libMesh::out << "|U("
624  << system.time + system.deltat
625  << ")|= "
626  << system.calculate_norm(*system.solution, 0, H1)
627  << std::endl;
628
629  libMesh::out << "|U("
630  << system.time
631  << ")|= "
632  << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
633  << std::endl;
634
635  libMesh::out << "|Z("
636  << system.time
637  << ")|= "
638  << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
639  << std::endl
640  << std::endl;
641
642  // Call the postprocess function which we have overloaded to compute
643  // accumulate the perturbed residuals
644  dynamic_cast<HeatSystem &>(system).perturb_accumulate_residuals(dynamic_cast<HeatSystem &>(system).get_parameter_vector());
645
646  // Now that we computed the accumulated, perturbed residuals, we can compute the
647  // approximate sensitivity
648  Number sensitivity_0_0 = (dynamic_cast<HeatSystem &>(system)).compute_final_sensitivity();
649
650  // Print it out
651  libMesh::out << "Sensitivity of QoI 0 w.r.t parameter 0 is: "
652  << sensitivity_0_0
653  << std::endl;
654
655  // Hard coded assert to ensure that the actual numbers we are
656  // getting are what they should be
657  // The 2e-4 tolerance is chosen to ensure success even with
658  // 32-bit floats
659  libmesh_assert_less(std::abs(sensitivity_0_0 - (-5.37173)), 2.e-4);
660
661 #ifdef NDEBUG
662  }
663  catch (...)
664  {
665  libMesh::err << '[' << mesh.processor_id()
666  << "] Caught exception; exiting early." << std::endl;
667  }
668 #endif
669
670  libMesh::err << '[' << mesh.processor_id()
671  << "] Completing output."
672  << std::endl;
673
674  // All done.
675  return 0;
676
677 #endif // LIBMESH_ENABLE_AMR
678 }
OStreamProxy err
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1545
bool print_solution_norms
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
void write_output(EquationSystems &es, unsigned int t_step, std::string solution_type, FEMParameters &param)
bool print_element_residuals
libMesh::Real initial_linear_tolerance
Real verify_analytic_jacobians
If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calc...
Definition: fem_system.h:203
void finish_initialization()
Definition: initial.C:11
Definition: system.h:378
bool analytic_jacobians
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Definition: diff_system.h:329
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
A member int that can be employed to indicate increased or reduced quadrature order.
Definition: system.h:1508
Real numerical_jacobian_h
If calculating numeric jacobians is required, the FEMSystem will perturb each solution vector entry b...
Definition: fem_system.h:175
std::string timesolver_core
Definition: femparameters.h:37
ElemType
Defines an enum for geometric element types.
bool require_residual_reduction
MeshBase & mesh
bool print_residual_norms
int main(int argc, char **argv)
libMesh::Real relative_step_tolerance
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr) const
Projects arbitrary functions onto the current solution.
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:334
std::vector< unsigned int > fe_order
bool print_element_residuals
Set print_element_residuals to true to print each R_elem contribution.
Definition: diff_system.h:344
Number adjoint_initial_value(const Point &p, const Parameters &, const std::string &, const std::string &)
PetscDiffSolver & solver
SolverPackage default_solver_package()
Definition: libmesh.C:995
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
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
const DofMap & get_dof_map() const
Definition: system.h:2030
unsigned int max_linear_iterations
int8_t boundary_id_type
Definition: id_types.h:51
bool print_residual_norms
Set print_residual_norms to true to print |F| whenever it is assembled.
Definition: diff_system.h:319
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
libMesh::Real minimum_linear_tolerance
std::vector< std::string > fe_family
libMesh::Real verify_analytic_jacobians
bool print_residuals
Set print_residuals to true to print F whenever it is assembled.
Definition: diff_system.h:324
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
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
Real calculate_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type, std::set< unsigned int > *skip_dimensions=libmesh_nullptr) const
Definition: system.C:1383
void libmesh_ignore(const T &)
void set_system_parameters(HeatSystem &system, FEMParameters &param)
void set_vector_preservation(const std::string &vec_name, bool preserve)
Allows one to set the boolean controlling whether the vector identified by vec_name should be "preser...
Definition: system.C:887
bool print_jacobian_norms
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
OStreamProxy out
libMesh::Real min_step_length
std::string & fe_family()
Definition: heatsystem.h:46
libMesh::Real timesolver_theta
Definition: femparameters.h:38
Gradient initial_grad(const Point &, const Parameters &, const std::string &, const std::string &)
Definition: initial.C:30
void read(GetPot &input, const std::vector< std::string > *other_variable_names=libmesh_nullptr)
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:989
bool print_element_jacobians
Set print_element_jacobians to true to print each J_elem contribution.
Definition: diff_system.h:349
bool print_element_jacobians
unsigned int max_nonlinear_iterations
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
std::map< libMesh::boundary_id_type, std::vector< unsigned int > > dirichlet_condition_variables
Definition: femparameters.h:90
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
Projects arbitrary functions onto a vector of degree of freedom values for the current system...
bool time_solver_quiet
bool & analytic_jacobians()
Definition: heatsystem.h:49