libMesh
heatsystem.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 
20 #include "libmesh/getpot.h"
21 
22 #include "heatsystem.h"
23 
24 #include "libmesh/fe_base.h"
25 #include "libmesh/fe_interface.h"
26 #include "libmesh/fem_context.h"
27 #include "libmesh/mesh.h"
28 #include "libmesh/quadrature.h"
29 #include "libmesh/string_to_enum.h"
30 #include "libmesh/system.h"
31 #include "libmesh/equation_systems.h"
32 #include "libmesh/zero_function.h"
33 #include "libmesh/dirichlet_boundaries.h"
34 #include "libmesh/dof_map.h"
35 
37 {
38  T_var = this->add_variable ("T", static_cast<Order>(_fe_order),
39  Utility::string_to_enum<FEFamily>(_fe_family));
40 
41  // Make sure the input file heat.in exists, and parse it.
42  {
43  std::ifstream i("heat.in");
44  if (!i)
45  libmesh_error_msg('[' << this->processor_id() << "] Can't find heat.in; exiting early.");
46  }
47  GetPot infile("heat.in");
48  _k = infile("k", 1.0);
49  _analytic_jacobians = infile("analytic_jacobians", true);
50 
51  parameters.push_back(_k);
52 
53  // Set equation system parameters _k and theta, so they can be read by the exact solution
54  this->get_equation_systems().parameters.set<Real>("_k") = _k;
55 
56  // The temperature is evolving, with a first order time derivative
57  this->time_evolving(T_var, 1);
58 
59  const boundary_id_type all_ids[4] = {0, 1, 2, 3};
60  std::set<boundary_id_type> all_bdys(all_ids, all_ids+(2*2));
61 
62  std::vector<unsigned int> T_only(1, T_var);
63 
64  ZeroFunction<Number> zero;
65 
66  // Most DirichletBoundary users will want to supply a "locally
67  // indexed" functor
69  (DirichletBoundary (all_bdys, T_only, zero,
71 
72  FEMSystem::init_data();
73 }
74 
75 
76 
77 void HeatSystem::init_context(DiffContext & context)
78 {
79  FEMContext & c = cast_ref<FEMContext &>(context);
80 
81  FEBase * elem_fe = libmesh_nullptr;
82  c.get_element_fe(0, elem_fe);
83 
84  // Now make sure we have requested all the data
85  // we need to build the linear system.
86  elem_fe->get_JxW();
87  elem_fe->get_dphi();
88 
89  // We'll have a more automatic solution to preparing adjoint
90  // solutions for time integration, eventually...
91  if (c.is_adjoint())
92  {
93  // A reference to the system context is built with
94  const System & sys = c.get_system();
95 
96  // Get a pointer to the adjoint solution vector
97  NumericVector<Number> & adjoint_solution =
98  const_cast<System &>(sys).get_adjoint_solution(0);
99 
100  // Add this adjoint solution to the vectors that diff context should localize
101  c.add_localized_vector(adjoint_solution, sys);
102  }
103 
104  FEMSystem::init_context(context);
105 }
106 
107 #define optassert(X) {if (!(X)) libmesh_error_msg("Assertion " #X " failed.");}
108 //#define optassert(X) libmesh_assert(X);
109 
110 bool HeatSystem::element_time_derivative (bool request_jacobian,
111  DiffContext & context)
112 {
113  bool compute_jacobian = request_jacobian && _analytic_jacobians;
114 
115  FEMContext & c = cast_ref<FEMContext &>(context);
116 
117  // First we get some references to cell-specific data that
118  // will be used to assemble the linear system.
119  FEBase * elem_fe = libmesh_nullptr;
120  c.get_element_fe(0, elem_fe);
121 
122  // Element Jacobian * quadrature weights for interior integration
123  const std::vector<Real> & JxW = elem_fe->get_JxW();
124 
125  // Element basis functions
126  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
127 
128  // Workaround for weird FC6 bug
129  optassert(c.get_dof_indices().size() > 0);
130 
131  // The number of local degrees of freedom in each variable
132  const unsigned int n_u_dofs = c.get_dof_indices(0).size();
133 
134  // The subvectors and submatrices we need to fill:
135  DenseSubMatrix<Number> & K = c.get_elem_jacobian(0, 0);
136  DenseSubVector<Number> & F = c.get_elem_residual(0);
137 
138  // Now we will build the element Jacobian and residual.
139  // Constructing the residual requires the solution and its
140  // gradient from the previous timestep. This must be
141  // calculated at each quadrature point by summing the
142  // solution degree-of-freedom values by the appropriate
143  // weight functions.
144  unsigned int n_qpoints = c.get_element_qrule().n_points();
145 
146  for (unsigned int qp=0; qp != n_qpoints; qp++)
147  {
148  // Compute the solution gradient at the Newton iterate
149  Gradient grad_T = c.interior_gradient(0, qp);
150 
151  for (unsigned int i=0; i != n_u_dofs; i++)
152  F(i) += JxW[qp] * -parameters[0] * (grad_T * dphi[i][qp]);
153  if (compute_jacobian)
154  for (unsigned int i=0; i != n_u_dofs; i++)
155  for (unsigned int j=0; j != n_u_dofs; ++j)
156  K(i,j) += JxW[qp] * -parameters[0] * (dphi[i][qp] * dphi[j][qp]);
157  } // end of the quadrature point qp-loop
158 
159  return compute_jacobian;
160 }
161 
162 // Perturb and accumulate dual weighted residuals
163 void HeatSystem::perturb_accumulate_residuals(ParameterVector & parameters_in)
164 {
165  const unsigned int Np = parameters_in.size();
166 
167  for (unsigned int j=0; j != Np; ++j)
168  {
169  Number old_parameter = *parameters_in[j];
170 
171  *parameters_in[j] = old_parameter - dp;
172 
173  this->assembly(true, false);
174 
175  this->rhs->close();
176 
177  UniquePtr<NumericVector<Number>> R_minus = this->rhs->clone();
178 
179  // The contribution at a single time step would be [f(z;p+dp) - <partialu/partialt, z>(p+dp) - <g(u),z>(p+dp)] * dt
180  // But since we compute the residual already scaled by dt, there is no need for the * dt
181  R_minus_dp += -R_minus->dot(this->get_adjoint_solution(0));
182 
183  *parameters_in[j] = old_parameter + dp;
184 
185  this->assembly(true, false);
186 
187  this->rhs->close();
188 
189  UniquePtr<NumericVector<Number>> R_plus = this->rhs->clone();
190 
191  R_plus_dp += -R_plus->dot(this->get_adjoint_solution(0));
192 
193  *parameters_in[j] = old_parameter;
194  }
195 }
Number R_plus_dp
Definition: heatsystem.h:132
Number R_minus_dp
Definition: heatsystem.h:133
ImplicitSystem & sys
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
Definition: heatsystem.C:110
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
const class libmesh_nullptr_t libmesh_nullptr
std::vector< Number > parameters
Definition: heatsystem.h:106
NumericVector< Number > * rhs
The system matrix.
const Number zero
.
Definition: libmesh.h:178
unsigned int T_var
Definition: heatsystem.h:126
const DofMap & get_dof_map() const
Definition: system.h:2030
int8_t boundary_id_type
Definition: id_types.h:51
bool _analytic_jacobians
Definition: heatsystem.h:129
virtual void init_context(DiffContext &context)
Definition: heatsystem.C:77
NumberVectorValue Gradient
virtual void time_evolving(unsigned int var)
Tells the DiffSystem that variable var is evolving with respect to time.
Definition: diff_physics.h:248
FEGenericBase< Real > FEBase
void perturb_accumulate_residuals(ParameterVector &parameters)
Definition: heatsystem.C:163
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Definition: parameters.h:469
std::string _fe_family
Definition: heatsystem.h:122
unsigned int _fe_order
Definition: heatsystem.h:123
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override
Prepares matrix or rhs for matrix assembly.
Definition: fem_system.C:852
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: heatsystem.C:36
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:989
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
System(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
Definition: system.C:59
processor_id_type processor_id() const
virtual UniquePtr< NumericVector< T > > clone() const =0