1 #include <libmesh/dof_map.h> 2 #include <libmesh/enum_solver_type.h> 3 #include <libmesh/enum_preconditioner_type.h> 4 #include <libmesh/fem_system.h> 5 #include <libmesh/newton_solver.h> 6 #include <libmesh/numeric_vector.h> 7 #include <libmesh/parallel.h> 16 template<
typename TimeSolverType>
27 template<
typename SystemType>
33 SystemType & system = es.
add_system<SystemType>(
"ScalarSystem");
35 system.time_solver = std::make_unique<TimeSolverType>(system);
39 DiffSolver & solver = *(system.time_solver->diff_solver().get());
44 NewtonSolver & newton = cast_ref<NewtonSolver &>(solver);
51 system.deltat = deltat;
53 TimeSolverType * time_solver = cast_ptr<TimeSolverType *>(system.time_solver.get());
54 this->aux_time_solver_init(*time_solver);
60 std::vector<dof_id_type> solution_index;
61 solution_index.push_back(0);
62 const bool has_solution = system.get_dof_map().all_semilocal_indices(solution_index);
64 for (
unsigned int t_step=0; t_step != n_timesteps; ++t_step)
67 system.time_solver->advance_timestep();
73 Number exact_soln = system.u(system.time);
74 rel_error =
std::abs((exact_soln - (*system.solution)(0))/exact_soln);
76 system.comm().max(rel_error);
79 LIBMESH_ASSERT_FP_EQUAL( rel_error,
81 std::numeric_limits<Real>::epsilon()*10 );
96 const std::string & name_in,
97 const unsigned int number_in)
114 this->time_evolving(_u_var,1);
122 FEMContext & c = cast_ref<FEMContext &>(context);
125 const unsigned int n_u_dofs =
129 for (
unsigned int qp=0; qp != n_qpoints; qp++)
131 Number Fval = this->F(c,qp);
133 for (
unsigned int i=0; i != n_u_dofs; i++)
137 return request_jacobian;
144 FEMContext & c = cast_ref<FEMContext &>(context);
148 const unsigned int n_u_dofs =
152 for (
unsigned int qp=0; qp != n_qpoints; qp++)
157 Number Mval = this->M(c,qp);
159 for (
unsigned int i=0; i != n_u_dofs; i++)
163 if (request_jacobian)
164 for (
unsigned int j=0; j != n_u_dofs; j++)
169 return request_jacobian;
186 const std::string & name_in,
187 const unsigned int number_in)
194 this->time_evolving(_u_var,2);
205 FEMContext & c = cast_ref<FEMContext &>(context);
209 const unsigned int n_u_dofs =
213 for (
unsigned int qp=0; qp != n_qpoints; qp++)
218 Number Cval = this->C(c,qp);
220 for (
unsigned int i=0; i != n_u_dofs; i++)
224 if (request_jacobian)
225 for (
unsigned int j=0; j != n_u_dofs; j++)
230 return request_jacobian;
236 FEMContext & c = cast_ref<FEMContext &>(context);
240 const unsigned int n_u_dofs =
244 for (
unsigned int qp=0; qp != n_qpoints; qp++)
249 Number Mval = this->M(c,qp);
251 for (
unsigned int i=0; i != n_u_dofs; i++)
255 if (request_jacobian)
256 for (
unsigned int j=0; j != n_u_dofs; j++)
261 return request_jacobian;
276 const std::string & name_in,
277 const unsigned int number_in)
285 FEMContext & c = cast_ref<FEMContext &>(context);
287 unsigned int v_var = this->get_second_order_dot_var(_u_var);
291 const unsigned int n_u_dofs =
295 for (
unsigned int qp=0; qp != n_qpoints; qp++)
297 Number Fval = this->F(c,qp);
299 for (
unsigned int i=0; i != n_u_dofs; i++)
303 return request_jacobian;
310 FEMContext & c = cast_ref<FEMContext &>(context);
312 unsigned int v_var = this->get_second_order_dot_var(_u_var);
318 const unsigned int n_u_dofs =
322 for (
unsigned int qp=0; qp != n_qpoints; qp++)
327 Number Cval = this->C(c,qp);
329 for (
unsigned int i=0; i != n_u_dofs; i++)
333 if (request_jacobian)
334 for (
unsigned int j=0; j != n_u_dofs; j++)
339 return request_jacobian;
345 FEMContext & c = cast_ref<FEMContext &>(context);
347 unsigned int v_var = this->get_second_order_dot_var(_u_var);
352 const unsigned int n_u_dofs =
356 for (
unsigned int qp=0; qp != n_qpoints; qp++)
361 Number Mval = this->M(c,qp);
363 for (
unsigned int i=0; i != n_u_dofs; i++)
367 if (request_jacobian)
368 for (
unsigned int j=0; j != n_u_dofs; j++)
373 return request_jacobian;
void run_test_with_exact_soln(Real deltat, unsigned int n_timesteps)
LinearSolver< Number > & get_linear_solver()
This is the EquationSystems class.
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
This class provides all data required for a physics package (e.g.
libMesh::Parallel::Communicator * TestCommWorld
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context) override
Note the nonlinear residual is M(u)
void interior_accel(unsigned int var, unsigned int qp, OutputType &u) const
Real absolute_residual_tolerance
The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolera...
SecondOrderScalarSystemSecondOrderTimeSolverBase(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
The libMesh namespace provides an interface to certain functionality in the library.
This class provides a specific system class.
virtual bool mass_residual(bool request_jacobian, DiffContext &context) override
Note the nonlinear residual is F(u)-M(u)
This class defines a solver which uses the default libMesh linear solver in a quasiNewton method to h...
Defines a dense subvector for use in finite element computations.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context) override
Note the nonlinear residual is F(u)-M(u)
Defines a dense submatrix for use in Finite Element-type computations.
void set_preconditioner_type(const PreconditionerType pct)
Sets the type of preconditioner to use.
virtual void init_data() override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
This class provides all data required for a physics package (e.g.
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
unsigned int n_points() const
void set_solver_type(const SolverType st)
Sets the type of solver to use.
FEMSystem-based class for testing of TimeSolvers using first order SCALARs.
Real get_elem_solution_accel_derivative() const
The derivative of the current elem_solution_accel w.r.t.
virtual bool mass_residual(bool request_jacobian, DiffContext &context) override
Note the nonlinear residual is F(u)-M(u)
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
FirstOrderScalarSystemBase(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
FEMSystem-based class for testing of TimeSolvers using second order SCALARs.
Real get_elem_solution_rate_derivative() const
The derivative of the current elem_solution_rate w.r.t.
SecondOrderScalarSystemFirstOrderTimeSolverBase(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
virtual void init()
Initialize all the systems.
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
virtual void init_data() override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
virtual bool mass_residual(bool request_jacobian, DiffContext &context) override
Note the nonlinear residual is F(u)-M(u)
void interior_rate(unsigned int var, unsigned int qp, OutputType &u) const
FEMSystem-based class for testing of TimeSolvers using second order SCALARs.
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
virtual bool damping_residual(bool request_jacobian, DiffContext &context) override
Note the nonlinear residual is M(u)
Real relative_residual_tolerance
Real relative_step_tolerance
virtual void init_data() override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
virtual void aux_time_solver_init(TimeSolverType &)
virtual bool damping_residual(bool request_jacobian, DiffContext &context) override
Note the nonlinear residual is M(u){u} + C(u)