libMesh
time_solver_test_common.h
Go to the documentation of this file.
1 #include <libmesh/dof_map.h>
2 #include <libmesh/fem_system.h>
3 
4 using namespace libMesh;
5 
6 template<typename TimeSolverType>
8 {
9 protected:
10 
11  // Any specialized initialization that's needed for the test
12  virtual void aux_time_solver_init( TimeSolverType & /*time_solver*/ ){}
13 
14  // Implementation for solving ODE of SystemType
15  // Note this test assumes that the time integrator gets the *exact* solution
16  // to within floating point tolerance.
17  template<typename SystemType>
18  void run_test_with_exact_soln(Real deltat, unsigned int n_timesteps)
19  {
22  EquationSystems es(mesh);
23  SystemType & system = es.add_system<SystemType>("ScalarSystem");
24 
25  system.time_solver.reset(new TimeSolverType(system));
26 
27  es.init();
28 
29  DiffSolver & solver = *(system.time_solver->diff_solver().get());
30  solver.relative_step_tolerance = std::numeric_limits<Real>::epsilon()*10;
31  solver.relative_residual_tolerance = std::numeric_limits<Real>::epsilon()*10;
32  solver.absolute_residual_tolerance = std::numeric_limits<Real>::epsilon()*10;
33 
34  system.deltat = deltat;
35 
36  TimeSolverType * time_solver = cast_ptr<TimeSolverType *>(system.time_solver.get());
37  this->aux_time_solver_init(*time_solver);
38 
39  // We're going to want to check our solution, and when we run
40  // "make check" with LIBMESH_RUN='mpirun -np N" for N>1 then we'll
41  // need to avoid checking on the processors that are just
42  // twiddling their thumbs, not owning our mesh point.
43  std::vector<dof_id_type> solution_index;
44  solution_index.push_back(0);
45  const bool has_solution = system.get_dof_map().all_semilocal_indices(solution_index);
46 
47  for (unsigned int t_step=0; t_step != n_timesteps; ++t_step)
48  {
49  system.solve();
50  system.time_solver->advance_timestep();
51 
52  if (has_solution)
53  {
54  // Use relative error for comparison, so "exact" is 0
55  Number exact_soln = system.u(system.time);
56  Real rel_error = std::abs((exact_soln - (*system.solution)(0))/exact_soln);
57 
58  CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0,
59  rel_error,
60  std::numeric_limits<Real>::epsilon()*10 );
61  }
62  }
63  }
64 
65 };
66 
68 
73 {
74 public:
76  const std::string & name_in,
77  const unsigned int number_in)
78  : FEMSystem(es, name_in, number_in)
79  {}
80 
81 
83  virtual Number F( FEMContext & context, unsigned int qp ) =0;
84 
86  virtual Number M( FEMContext & context, unsigned int qp ) =0;
87 
89  virtual Number u( Real t ) =0;
90 
91  virtual void init_data () libmesh_override
92  {
93  _u_var = this->add_variable ("u", FIRST, LAGRANGE);
94  this->time_evolving(_u_var,1);
96  }
97 
99  virtual bool element_time_derivative (bool request_jacobian,
100  DiffContext & context) libmesh_override
101  {
102  FEMContext & c = cast_ref<FEMContext &>(context);
103  DenseSubVector<Number> & Fu = c.get_elem_residual(_u_var);
104 
105  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
106  unsigned int n_qpoints = c.get_element_qrule().n_points();
107 
108  for (unsigned int qp=0; qp != n_qpoints; qp++)
109  {
110  Number Fval = this->F(c,qp);
111 
112  for (unsigned int i=0; i != n_u_dofs; i++)
113  Fu(i) += Fval;
114  }
115 
116  return request_jacobian;
117  }
118 
120  virtual bool mass_residual (bool request_jacobian,
121  DiffContext & context) libmesh_override
122  {
123  FEMContext & c = cast_ref<FEMContext &>(context);
124  DenseSubVector<Number> & Fu = c.get_elem_residual(_u_var);
125  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(_u_var, _u_var);
126 
127  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
128  unsigned int n_qpoints = c.get_element_qrule().n_points();
129 
130  for (unsigned int qp=0; qp != n_qpoints; qp++)
131  {
132  libMesh::Number udot;
133  c.interior_rate(_u_var,qp,udot);
134 
135  Number Mval = this->M(c,qp);
136 
137  for (unsigned int i=0; i != n_u_dofs; i++)
138  {
139  Fu(i) -= Mval*udot;
140 
141  if (request_jacobian)
142  for (unsigned int j=0; j != n_u_dofs; j++)
143  Kuu(i,j) -= Mval*context.get_elem_solution_rate_derivative();
144  }
145  }
146 
147  return request_jacobian;
148  }
149 
150 protected:
151  unsigned int _u_var;
152 };
153 
155 
161 {
162 public:
164  const std::string & name_in,
165  const unsigned int number_in)
166  : FirstOrderScalarSystemBase(es, name_in, number_in)
167  {}
168 
169  virtual void init_data () libmesh_override
170  {
171  _u_var = this->add_variable ("u", FIRST, LAGRANGE);
172  this->time_evolving(_u_var,2);
174  }
175 
177  virtual Number C( FEMContext & context, unsigned int qp ) =0;
178 
180  virtual bool damping_residual (bool request_jacobian,
181  DiffContext & context) libmesh_override
182  {
183  FEMContext & c = cast_ref<FEMContext &>(context);
184  DenseSubVector<Number> & Fu = c.get_elem_residual(_u_var);
185  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(_u_var, _u_var);
186 
187  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
188  unsigned int n_qpoints = c.get_element_qrule().n_points();
189 
190  for (unsigned int qp=0; qp != n_qpoints; qp++)
191  {
192  libMesh::Number udot;
193  c.interior_rate(_u_var,qp,udot);
194 
195  Number Cval = this->C(c,qp);
196 
197  for (unsigned int i=0; i != n_u_dofs; i++)
198  {
199  Fu(i) += Cval*udot;
200 
201  if (request_jacobian)
202  for (unsigned int j=0; j != n_u_dofs; j++)
203  Kuu(i,j) += Cval*context.get_elem_solution_rate_derivative();
204  }
205  }
206 
207  return request_jacobian;
208  }
209 
210  virtual bool mass_residual (bool request_jacobian,
211  DiffContext & context) libmesh_override
212  {
213  FEMContext & c = cast_ref<FEMContext &>(context);
214  DenseSubVector<Number> & Fu = c.get_elem_residual(_u_var);
215  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(_u_var, _u_var);
216 
217  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
218  unsigned int n_qpoints = c.get_element_qrule().n_points();
219 
220  for (unsigned int qp=0; qp != n_qpoints; qp++)
221  {
222  libMesh::Number uddot;
223  c.interior_accel(_u_var,qp,uddot);
224 
225  Number Mval = this->M(c,qp);
226 
227  for (unsigned int i=0; i != n_u_dofs; i++)
228  {
229  Fu(i) += Mval*uddot;
230 
231  if (request_jacobian)
232  for (unsigned int j=0; j != n_u_dofs; j++)
233  Kuu(i,j) += Mval*context.get_elem_solution_accel_derivative();
234  }
235  }
236 
237  return request_jacobian;
238  }
239 };
240 
241 
243 
249 {
250 public:
252  const std::string & name_in,
253  const unsigned int number_in)
254  : SecondOrderScalarSystemSecondOrderTimeSolverBase(es, name_in, number_in)
255  {}
256 
258  virtual bool element_time_derivative (bool request_jacobian,
259  DiffContext & context) libmesh_override
260  {
261  FEMContext & c = cast_ref<FEMContext &>(context);
262 
263  unsigned int v_var = this->get_second_order_dot_var(_u_var);
264 
266 
267  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
268  unsigned int n_qpoints = c.get_element_qrule().n_points();
269 
270  for (unsigned int qp=0; qp != n_qpoints; qp++)
271  {
272  Number Fval = this->F(c,qp);
273 
274  for (unsigned int i=0; i != n_u_dofs; i++)
275  Fv(i) += Fval;
276  }
277 
278  return request_jacobian;
279  }
280 
282  virtual bool damping_residual (bool request_jacobian,
283  DiffContext & context) libmesh_override
284  {
285  FEMContext & c = cast_ref<FEMContext &>(context);
286 
287  unsigned int v_var = this->get_second_order_dot_var(_u_var);
288 
290 
291  DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
292 
293  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
294  unsigned int n_qpoints = c.get_element_qrule().n_points();
295 
296  for (unsigned int qp=0; qp != n_qpoints; qp++)
297  {
298  libMesh::Number udot;
299  c.interior_rate(v_var,qp,udot);
300 
301  Number Cval = this->C(c,qp);
302 
303  for (unsigned int i=0; i != n_u_dofs; i++)
304  {
305  Fv(i) += Cval*udot;
306 
307  if (request_jacobian)
308  for (unsigned int j=0; j != n_u_dofs; j++)
309  Kvv(i,j) += Cval*context.get_elem_solution_rate_derivative();
310  }
311  }
312 
313  return request_jacobian;
314  }
315 
316  virtual bool mass_residual (bool request_jacobian,
317  DiffContext & context) libmesh_override
318  {
319  FEMContext & c = cast_ref<FEMContext &>(context);
320 
321  unsigned int v_var = this->get_second_order_dot_var(_u_var);
322 
324  DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
325 
326  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
327  unsigned int n_qpoints = c.get_element_qrule().n_points();
328 
329  for (unsigned int qp=0; qp != n_qpoints; qp++)
330  {
331  libMesh::Number uddot;
332  c.interior_accel(v_var,qp,uddot);
333 
334  Number Mval = this->M(c,qp);
335 
336  for (unsigned int i=0; i != n_u_dofs; i++)
337  {
338  Fv(i) += Mval*uddot;
339 
340  if (request_jacobian)
341  for (unsigned int j=0; j != n_u_dofs; j++)
342  Kvv(i,j) += Mval*context.get_elem_solution_accel_derivative();
343  }
344  }
345 
346  return request_jacobian;
347  }
348 };
void run_test_with_exact_soln(Real deltat, unsigned int n_timesteps)
double abs(double a)
This is the EquationSystems class.
void build_point(UnstructuredMesh &mesh, const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 0D meshes.
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:28
virtual bool mass_residual(bool request_jacobian, DiffContext &context) libmesh_override
Note the nonlinear residual is F(u)-M(u)
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
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
SecondOrderScalarSystemSecondOrderTimeSolverBase(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
MeshBase & mesh
virtual bool damping_residual(bool request_jacobian, DiffContext &context) libmesh_override
Note the nonlinear residual is M(u)
virtual bool mass_residual(bool request_jacobian, DiffContext &context) libmesh_override
Note the nonlinear residual is F(u)-M(u)
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context) libmesh_override
Note the nonlinear residual is M(u)
The libMesh namespace provides an interface to certain functionality in the library.
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
virtual bool mass_residual(bool request_jacobian, DiffContext &context) libmesh_override
Note the nonlinear residual is F(u)-M(u)
This class provides a specific system class.
Definition: fem_system.h:53
Real get_elem_solution_accel_derivative() const
The derivative of the current elem_solution_accel w.r.t.
Definition: diff_context.h:437
void interior_rate(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1273
Defines a dense subvector for use in finite element computations.
PetscDiffSolver & solver
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
Definition: diff_solver.h:70
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
Real get_elem_solution_rate_derivative() const
The derivative of the current elem_solution_rate w.r.t.
Definition: diff_context.h:428
Defines a dense submatrix for use in Finite Element-type computations.
void interior_accel(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1295
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
FEMSystem-based class for testing of TimeSolvers using first order SCALARs.
FirstOrderScalarSystemBase(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool damping_residual(bool request_jacobian, DiffContext &context) libmesh_override
Note the nonlinear residual is M(u){u} + C(u)
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
FEMSystem-based class for testing of TimeSolvers using second order SCALARs.
SecondOrderScalarSystemFirstOrderTimeSolverBase(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
virtual void init()
Initialize all the systems.
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:845
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
FEMSystem-based class for testing of TimeSolvers using second order SCALARs.
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context) libmesh_override
Note the nonlinear residual is F(u)-M(u)
unsigned int n_points() const
Definition: quadrature.h:113
Real relative_residual_tolerance
Definition: diff_solver.h:192
virtual void aux_time_solver_init(TimeSolverType &)