Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 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_context.h"
      20             : #include "libmesh/diff_physics.h"
      21             : #include "libmesh/system.h"
      22             : 
      23             : namespace libMesh
      24             : {
      25             : 
      26       71800 : DifferentiablePhysics::~DifferentiablePhysics() = default;
      27             : 
      28           0 : void DifferentiablePhysics::clear_physics ()
      29             : {
      30           0 :   _time_evolving.resize(0);
      31           0 : }
      32             : 
      33             : 
      34             : 
      35        5044 : void DifferentiablePhysics::init_physics (const System & sys)
      36             : {
      37             :   // give us flags for every variable that might be time evolving
      38        5194 :   _time_evolving.resize(sys.n_vars(), false);
      39        5044 : }
      40             : 
      41        7166 : void DifferentiablePhysics::time_evolving (unsigned int var,
      42             :                                            unsigned int order)
      43             : {
      44        7166 :   libmesh_error_msg_if(order != 1 && order != 2, "Input order must be 1 or 2!");
      45             : 
      46        7378 :   if (_time_evolving.size() <= var)
      47        6430 :     _time_evolving.resize(var+1, 0);
      48             : 
      49        7166 :   _time_evolving[var] = order;
      50             : 
      51        7166 :   if (order == 1)
      52        4884 :     _first_order_vars.insert(var);
      53             :   else
      54        2070 :     _second_order_vars.insert(var);
      55        7166 : }
      56             : 
      57           0 : bool DifferentiablePhysics::nonlocal_mass_residual(bool request_jacobian,
      58             :                                                    DiffContext & c)
      59             : {
      60           0 :   FEMContext & context = cast_ref<FEMContext &>(c);
      61             : 
      62           0 :   for (auto var : make_range(context.n_vars()))
      63             :     {
      64           0 :       if (!this->is_time_evolving(var))
      65           0 :         continue;
      66             : 
      67           0 :       if (c.get_system().variable(var).type().family != SCALAR)
      68           0 :         continue;
      69             : 
      70             :       const std::vector<dof_id_type> & dof_indices =
      71           0 :         context.get_dof_indices(var);
      72             : 
      73             :       const unsigned int n_dofs = cast_int<unsigned int>
      74           0 :         (dof_indices.size());
      75             : 
      76           0 :       DenseSubVector<Number> & Fs = context.get_elem_residual(var);
      77           0 :       DenseSubMatrix<Number> & Kss = context.get_elem_jacobian( var, var );
      78             : 
      79             :       const libMesh::DenseSubVector<libMesh::Number> & Us =
      80           0 :         context.get_elem_solution(var);
      81             : 
      82           0 :       for (unsigned int i=0; i != n_dofs; ++i)
      83             :         {
      84           0 :           Fs(i) -= Us(i);
      85             : 
      86           0 :           if (request_jacobian)
      87           0 :             Kss(i,i) -= context.elem_solution_rate_derivative;
      88             :         }
      89             :     }
      90             : 
      91           0 :   return request_jacobian;
      92             : }
      93             : 
      94             : 
      95             : 
      96    20563472 : bool DifferentiablePhysics::_eulerian_time_deriv (bool request_jacobian,
      97             :                                                   DiffContext & context)
      98             : {
      99             :   // For any problem we need time derivative terms
     100             :   request_jacobian =
     101    20563472 :     this->element_time_derivative(request_jacobian, context);
     102             : 
     103             :   // For a moving mesh problem we may need the pseudoconvection term too
     104    20563472 :   return this->eulerian_residual(request_jacobian, context) &&
     105    20563472 :     request_jacobian;
     106             : }
     107             : 
     108             : } // namespace libMesh