1 // The libMesh Finite Element Library. 2 // Copyright (C) 2002-2020 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 #ifndef LIBMESH_FIRST_ORDER_UNSTEADY_SOLVER_H 19 #define LIBMESH_FIRST_ORDER_UNSTEADY_SOLVER_H 20 21 #include "libmesh/unsteady_solver.h" 22 23 namespace libMesh 24 { 25 /** 26 * Generic class from which first order UnsteadySolvers should subclass. 27 * 28 * Subclasses of this class are meant to solve problems of the form 29 * \f[ M(u)\dot{u} = F(u)\f] 30 * 31 * There is also infrastructure for subclasses to support solving second 32 * order-in-time systems (or both first order and second order). 33 * In particular, consider the second system 34 * \f[ M(u)\ddot{u} + C(u)\dot{u} + F(u) = 0 \f] 35 * If we introduce the equation \f$\dot{u} = v\f$, then, we have the 36 * following first order system: 37 * \f{eqnarray}{ \dot{u} &=& v \\ M(u)\dot{v} + C(u)v + F(u) &=& 0\f} 38 * Such systems are supported by the user specifying that the time order 39 * of the variable \f$u\f$ is 2 when adding the variable using 40 * FEMSystem::time_evolving. This class will then add the \f$v\f$ variables 41 * to the FEMSystem (named "dot_u" if the second order variable name is "u"). 42 * Subclasses will then need to make sure to use the function 43 * FirstOrderUnsteadySolver::prepare_accel() to prepare data structures for the 44 * users to use from FEMContext. Furthermore, subclasses will need to appropriately 45 * call damping_residual functions. Finally, subclasses will call 46 * FirstOrderUnsteadySolver::compute_second_order_eqns() to actually assemble 47 * residual and Jacobians for the \f$ \dot{u} = v\f$ equations. These aspects 48 * should be invisible ot users during their element assembly. 49 * 50 * Unfortunately, complete usage of the feature of treating second order 51 * equations as a system of first order equations is not completely invisible 52 * to the user. The user must assemble their residual in the correct equation. 53 * In particular, the residual must go in the \f$v\f$ residual equation, i.e. 54 * \f$ M(u)\dot{v} + C(u)v + F(u) = 0 \f$, and, 55 * subsequently, they must also be careful to populate to the correct Jacobian 56 * blocks. The function to help facilitate this is get_second_order_dot_var. 57 * If you have a second order variable, you pass that variable index and 58 * get_second_order_dot_var will return the index of the corresponding velocity 59 * variable, \f$v\f$ in the notation above. Then, the user knows what blocks to 60 * populate. 61 * 62 * \note The API is designed so that if the user codes to this 63 * paradigm, the TimeSolver will be interchangeable for those element 64 * kernels. That is, they'll be able to use either a 65 * FirstOrderUnsteadySolver or a SecondOrderUnsteadySolver. 66 * 67 * This class is part of the new DifferentiableSystem framework, 68 * which is still experimental. Users of this framework should 69 * beware of bugs and future API changes. 70 * 71 * \author Paul T. Bauman 72 * \date 2015 73 */ 74 class FirstOrderUnsteadySolver : public UnsteadySolver 75 { 76 public: 77 /** 78 * Constructor. Requires a reference to the system 79 * to be solved. 80 */ 81 explicit FirstOrderUnsteadySolver(sys_type & s)82 FirstOrderUnsteadySolver (sys_type & s) 83 : UnsteadySolver(s) {} 84 85 /** 86 * Destructor. 87 */ ~FirstOrderUnsteadySolver()88 virtual ~FirstOrderUnsteadySolver (){} 89 time_order()90 virtual unsigned int time_order() const override 91 { return 1; } 92 93 protected: 94 95 /** 96 * If there are second order variables in the system, 97 * then we also prepare the accel for those variables 98 * so the user can treat them as such. 99 */ 100 void prepare_accel(DiffContext & context); 101 102 /** 103 * If there are second order variables, then we need to compute their residual equations 104 * and corresponding Jacobian. The residual equation will simply be 105 * \f$ \dot{u} - v = 0 \f$, where \f$ u \f$ is the second order variable add 106 * by the user and \f$ v \f$ is the variable added by the time-solver as the 107 * "velocity" variable. 108 */ 109 bool compute_second_order_eqns(bool compute_jacobian, DiffContext & c); 110 111 }; 112 113 } // end namespace libMesh 114 115 #endif // LIBMESH_FIRST_ORDER_UNSTEADY_SOLVER_H 116