libMesh
diff_system.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 #include "libmesh/diff_solver.h"
20 #include "libmesh/diff_system.h"
21 #include "libmesh/time_solver.h"
22 #include "libmesh/unsteady_solver.h"
23 #include "libmesh/dirichlet_boundaries.h"
24 #include "libmesh/dof_map.h"
25 #include "libmesh/zero_function.h"
26 
27 #include <utility> // std::swap
28 
29 namespace libMesh
30 {
31 
32 
33 
35  const std::string & name_in,
36  const unsigned int number_in) :
37  Parent (es, name_in, number_in),
38  time_solver (),
39  deltat(1.),
40  print_solution_norms(false),
41  print_solutions(false),
42  print_residual_norms(false),
43  print_residuals(false),
44  print_jacobian_norms(false),
45  print_jacobians(false),
46  print_element_solutions(false),
47  print_element_residuals(false),
48  print_element_jacobians(false),
49  _diff_physics(this),
50  diff_qoi(this)
51 {
52 }
53 
54 
55 
57 {
58  // If we had an attached Physics object, delete it.
59  if (this->_diff_physics != this)
60  delete this->_diff_physics;
61 
62  // If we had an attached QoI object, delete it.
63  if (this->diff_qoi != this)
64  delete this->diff_qoi;
65 }
66 
67 
68 
70 {
71  // If we had an attached Physics object, delete it.
72  if (this->_diff_physics != this)
73  {
74  delete this->_diff_physics;
75  this->_diff_physics = this;
76  }
77  // If we had no attached Physics object, clear our own Physics data
78  else
79  this->clear_physics();
80 
81  // If we had an attached QoI object, delete it.
82  if (this->diff_qoi != this)
83  {
84  delete this->diff_qoi;
85  this->diff_qoi = this;
86  }
87  // If we had no attached QoI object, clear our own QoI data
88  else
89  this->clear_qoi();
90 
91  use_fixed_solution = false;
92 }
93 
94 
95 
97 {
99 
101  libmesh_assert_equal_to (&(time_solver->system()), this);
102 
103  time_solver->reinit();
104 }
105 
106 
107 
109 {
110  // If it isn't a separate initialized-upon-attachment object, do any
111  // initialization our physics needs.
112  if (this->_diff_physics == this)
113  this->init_physics(*this);
114 
115  // Do any initialization our solvers need
117  libmesh_assert_equal_to (&(time_solver->system()), this);
118 
119  // Now check for second order variables and add their velocities to the System.
120  if (!time_solver->is_steady())
121  {
122  const UnsteadySolver & unsteady_solver =
123  cast_ref<const UnsteadySolver &>(*(time_solver.get()));
124 
125  if (unsteady_solver.time_order() == 1)
127  }
128 
129  time_solver->init();
130 
131  // Next initialize ImplicitSystem data
133 
134  time_solver->init_data();
135 }
136 
138 {
139  DiffContext * context = new DiffContext(*this);
140  context->set_deltat_pointer( &this->deltat );
141  return UniquePtr<DiffContext>(context);
142 }
143 
144 
146 {
147  this->assembly(true, true);
148 }
149 
150 
151 
153 {
154  // Get the time solver object associated with the system, and tell it that
155  // we are not solving the adjoint problem
156  this->get_time_solver().set_is_adjoint(false);
157 
158  libmesh_assert_equal_to (&(time_solver->system()), this);
159  time_solver->solve();
160 }
161 
162 
163 
164 std::pair<unsigned int, Real> DifferentiableSystem::adjoint_solve (const QoISet & qoi_indices)
165 {
166  // Get the time solver object associated with the system, and tell it that
167  // we are solving the adjoint problem
168  this->get_time_solver().set_is_adjoint(true);
169 
170  return this->ImplicitSystem::adjoint_solve(qoi_indices);
171 }
172 
173 
174 
176 {
178  libmesh_assert_equal_to (&(time_solver->system()), this);
179  return this->time_solver->linear_solver().get();
180 }
181 
182 
183 
184 std::pair<unsigned int, Real> DifferentiableSystem::get_linear_solve_parameters() const
185 {
187  libmesh_assert_equal_to (&(time_solver->system()), this);
188  return std::make_pair(this->time_solver->diff_solver()->max_linear_iterations,
189  this->time_solver->diff_solver()->relative_residual_tolerance);
190 }
191 
192 
193 
195 {
196 }
197 
199 {
200  const std::set<unsigned int> & second_order_vars = this->get_second_order_vars();
201  if (!second_order_vars.empty())
202  {
203  for (std::set<unsigned int>::const_iterator var_it = second_order_vars.begin();
204  var_it != second_order_vars.end(); ++var_it)
205  {
206  const Variable & var = this->variable(*var_it);
207  std::string new_var_name = std::string("dot_")+var.name();
208 
209  unsigned int v_var_idx;
210 
211  if (var.active_subdomains().empty())
212  v_var_idx = this->add_variable( new_var_name, var.type() );
213  else
214  v_var_idx = this->add_variable( new_var_name, var.type(), &var.active_subdomains() );
215 
216  _second_order_dot_vars.insert( std::pair<unsigned int,unsigned int>(*var_it,v_var_idx) );
217 
218  // The new velocities are time evolving variables of first order
219  this->time_evolving( v_var_idx, 1 );
220 
221  // And if there are any boundary conditions set on the second order
222  // variable, we also need to set it on its velocity variable.
223  this->add_dot_var_dirichlet_bcs( *var_it, v_var_idx );
224  }
225  }
226 }
227 
229  unsigned int dot_var_idx )
230 {
231  // We're assuming that there could be a lot more variables than
232  // boundary conditions, so we search each of the boundary conditions
233  // for this variable rather than looping over boundary conditions
234  // in a separate loop and searching through all the variables.
235  const DirichletBoundaries * all_dbcs =
237 
238  if (all_dbcs)
239  {
240  // We need to cache the DBCs to be added so that we add them
241  // after looping over the existing DBCs. Otherwise, we're polluting
242  // the thing we're looping over.
243  std::vector<DirichletBoundary*> new_dbcs;
244 
245  DirichletBoundaries::const_iterator dbc_it = all_dbcs->begin();
246  for ( ; dbc_it != all_dbcs->end(); ++dbc_it )
247  {
248  libmesh_assert(*dbc_it);
249  DirichletBoundary & dbc = *(*dbc_it);
250 
251  // Look for second order variable in the current
252  // DirichletBoundary object
253  std::vector<unsigned int>::const_iterator dbc_var_it =
254  std::find( dbc.variables.begin(), dbc.variables.end(), var_idx );
255 
256  // If we found it, then we also need to add it's corresponding
257  // "dot" variable to a DirichletBoundary
258  std::vector<unsigned int> vars_to_add;
259  if (dbc_var_it != dbc.variables.end())
260  vars_to_add.push_back(dot_var_idx);
261 
262  if (!vars_to_add.empty())
263  {
264  // We need to check if the boundary condition is time-dependent.
265  // Currently, we cannot automatically differentiate w.r.t. time
266  // so if the user supplies a time-dependent Dirichlet BC, then
267  // we can't automatically support the Dirichlet BC for the
268  // "velocity" boundary condition, so we error. Otherwise,
269  // the "velocity boundary condition will just be zero.
270  bool is_time_evolving_bc = false;
271  if (dbc.f)
272  is_time_evolving_bc = dbc.f->is_time_dependent();
273  else if (dbc.f_fem)
274  // We it's a FEMFunctionBase object, it will be implicitly
275  // time-dependent since it is assumed to depend on the solution.
276  is_time_evolving_bc = true;
277  else
278  libmesh_error_msg("Could not find valid boundary function!");
279 
280  if (is_time_evolving_bc)
281  libmesh_error_msg("Cannot currently support time-dependent Dirichlet BC for dot variables!");
282 
283 
284  DirichletBoundary * new_dbc;
285 
286  if (dbc.f)
287  {
289 
290  new_dbc = new DirichletBoundary(dbc.b, vars_to_add, zf);
291  }
292  else
293  libmesh_error();
294 
295  new_dbcs.push_back(new_dbc);
296  }
297  }
298 
299  // Now add the new DBCs for the "dot" vars to the DofMap
300  std::vector<DirichletBoundary*>::iterator new_dbc_it =
301  new_dbcs.begin();
302 
303  for ( ; new_dbc_it != new_dbcs.end(); ++new_dbc_it )
304  {
305  const DirichletBoundary & dbc = *(*new_dbc_it);
306  this->get_dof_map().add_dirichlet_boundary(dbc);
307  delete *new_dbc_it;
308  }
309 
310  } // if (all_dbcs)
311 }
312 
313 unsigned int DifferentiableSystem::get_second_order_dot_var( unsigned int var ) const
314 {
315  // For SteadySolver or SecondOrderUnsteadySolvers, we just give back var
316  unsigned int dot_var = var;
317 
318  if (!time_solver->is_steady())
319  {
320  const UnsteadySolver & unsteady_solver =
321  cast_ref<const UnsteadySolver &>(*(time_solver.get()));
322 
323  if (unsteady_solver.time_order() == 1)
324  dot_var = this->_second_order_dot_vars.find(var)->second;
325  }
326 
327  return dot_var;
328 }
329 
331 {
332  bool have_first_order_scalar_vars = false;
333 
334  if (this->have_first_order_vars())
335  {
336  for (std::set<unsigned int>::const_iterator var_it = this->get_first_order_vars().begin();
337  var_it != this->get_first_order_vars().end();
338  ++var_it)
339  {
340  if (this->variable(*var_it).type().family == SCALAR)
341  have_first_order_scalar_vars = true;
342  }
343  }
344 
346 }
347 
349 {
350  bool have_second_order_scalar_vars = false;
351 
352  if (this->have_second_order_vars())
353  {
354  for (std::set<unsigned int>::const_iterator var_it = this->get_second_order_vars().begin();
355  var_it != this->get_second_order_vars().end();
356  ++var_it)
357  {
358  if (this->variable(*var_it).type().family == SCALAR)
359  have_second_order_scalar_vars = true;
360  }
361  }
362 
364 }
365 
366 
367 
369 {
370  std::swap(this->_diff_physics, swap_physics);
371 }
372 
373 
374 } // namespace libMesh
FEFamily family
The type of finite element.
Definition: fe_type.h:203
const FEType & type() const
Definition: variable.h:119
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
This is the EquationSystems class.
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) libmesh_override
Assembles & solves the linear system (dR/du)^T*z = dq/du, for those quantities of interest q specifie...
const std::set< unsigned int > & get_second_order_vars() const
Definition: diff_physics.h:529
virtual void clear_physics()
Clear any data structures associated with the physics.
Definition: diff_physics.C:33
const std::string & name() const
Definition: variable.h:100
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
ConstFunction that simply returns 0.
Definition: zero_function.h:35
virtual void assemble() libmesh_override
Prepares matrix and rhs for matrix assembly.
Definition: diff_system.C:145
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
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
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: diff_system.C:108
virtual UniquePtr< DiffContext > build_context()
Builds a DiffContext object with enough information to do evaluations on each element.
Definition: diff_system.C:137
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
void add_second_order_dot_vars()
Helper function to add "velocity" variables that are cousins to second order-in-time variables in the...
Definition: diff_system.C:198
std::vector< unsigned int > variables
void add_dot_var_dirichlet_bcs(unsigned int var_idx, unsigned int dot_var_idx)
Helper function to and Dirichlet boundary conditions to "dot" variable cousins of second order variab...
Definition: diff_system.C:228
The libMesh namespace provides an interface to certain functionality in the library.
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
Definition: diff_system.C:175
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
Definition: diff_system.C:69
std::map< unsigned int, unsigned int > _second_order_dot_vars
If the user adds any second order variables, then we need to also cache the map to their correspondin...
Definition: diff_physics.h:570
bool have_second_order_scalar_vars() const
Check for any second order vars that are also belong to FEFamily::SCALAR.
Definition: diff_system.C:348
bool have_second_order_vars() const
Definition: diff_physics.h:523
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: diff_system.C:152
libmesh_assert(j)
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override=0
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
DifferentiableQoI * diff_qoi
Pointer to object to use for quantity of interest assembly evaluations.
Definition: diff_system.h:365
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual void init_physics(const System &sys)
Initialize any data structures associated with the physics.
Definition: diff_physics.C:40
unsigned int get_second_order_dot_var(unsigned int var) const
For a given second order (in time) variable var, this method will return the index to the correspondi...
Definition: diff_system.C:313
void swap_physics(DifferentiablePhysics *&swap_physics)
Swap current physics object with external object.
Definition: diff_system.C:368
DifferentiablePhysics * _diff_physics
Pointer to object to use for physics assembly evaluations.
Definition: diff_system.h:358
void set_is_adjoint(bool _is_adjoint_value)
Accessor for setting whether we need to do a primal or adjoint solve.
Definition: time_solver.h:239
This class defines the notion of a variable in the system.
Definition: variable.h:49
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual unsigned int time_order() const =0
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
Definition: system.h:2114
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
bool use_fixed_solution
A boolean to be set to true by systems using elem_fixed_solution, for optional use by e...
Definition: system.h:1493
virtual void release_linear_solver(LinearSolver< Number > *) const libmesh_override
Releases a pointer to a linear solver acquired by this->get_linear_solver()
Definition: diff_system.C:194
bool have_first_order_scalar_vars() const
Check for any first order vars that are also belong to FEFamily::SCALAR.
Definition: diff_system.C:330
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
const DirichletBoundaries * get_dirichlet_boundaries() const
Definition: dof_map.h:1145
virtual void clear_qoi()
Clear all the data structures associated with the QoI.
Definition: diff_qoi.h:74
virtual std::pair< unsigned int, Real > get_linear_solve_parameters() const libmesh_override
Definition: diff_system.C:184
virtual void time_evolving(unsigned int var)
Tells the DiffSystem that variable var is evolving with respect to time.
Definition: diff_physics.h:248
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
This class provides a specific system class.
Definition: diff_physics.h:74
DifferentiableSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
Definition: diff_system.C:34
UniquePtr< FunctionBase< Number > > f
This is a generic class that defines a solver to handle time integration of DifferentiableSystems.
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
const std::set< unsigned int > & get_first_order_vars() const
Definition: diff_physics.h:516
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
void set_deltat_pointer(Real *dt)
Points the _deltat member of this class at a timestep value stored in the creating System...
Definition: diff_context.C:138
const std::set< subdomain_id_type > & active_subdomains() const
Definition: variable.h:150
std::set< boundary_id_type > b
UniquePtr< FEMFunctionBase< Number > > f_fem
virtual ~DifferentiableSystem()
Destructor.
Definition: diff_system.C:56
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
Definition: diff_system.C:96
TimeSolver & get_time_solver()
Definition: diff_system.h:402
This class provides a specific system class.