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             : 
      20             : #include "libmesh/dof_map.h"
      21             : #include "libmesh/elem.h"
      22             : #include "libmesh/equation_systems.h"
      23             : #include "libmesh/fe_base.h"
      24             : #include "libmesh/fem_context.h"
      25             : #include "libmesh/fem_system.h"
      26             : #include "libmesh/libmesh_logging.h"
      27             : #include "libmesh/mesh_base.h"
      28             : #include "libmesh/numeric_vector.h"
      29             : #include "libmesh/quadrature.h"
      30             : #include "libmesh/sparse_matrix.h"
      31             : #include "libmesh/time_solver.h"
      32             : #include "libmesh/unsteady_solver.h" // For eulerian_residual
      33             : 
      34             : 
      35             : namespace libMesh
      36             : {
      37             : 
      38    20563472 : bool FEMPhysics::eulerian_residual (bool request_jacobian,
      39             :                                     DiffContext & /*c*/)
      40             : {
      41             :   // Only calculate a mesh movement residual if it's necessary
      42    20563472 :   if (!_mesh_sys)
      43     1845168 :     return request_jacobian;
      44             : 
      45           0 :   libmesh_not_implemented();
      46             : 
      47             : #if 0
      48             :   FEMContext & context = cast_ref<FEMContext &>(c);
      49             : 
      50             :   // This function only supports fully coupled mesh motion for now
      51             :   libmesh_assert_equal_to (_mesh_sys, this);
      52             : 
      53             :   unsigned int n_qpoints = (context.get_element_qrule())->n_points();
      54             : 
      55             :   const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ?
      56             :     0 : context.dof_indices_var[_mesh_x_var].size();
      57             :   const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ?
      58             :     0 : context.dof_indices_var[_mesh_y_var].size();
      59             :   const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ?
      60             :     0 : context.dof_indices_var[_mesh_z_var].size();
      61             : 
      62             :   const unsigned int mesh_xyz_var = n_x_dofs ? _mesh_x_var :
      63             :     (n_y_dofs ? _mesh_y_var :
      64             :      (n_z_dofs ? _mesh_z_var :
      65             :       libMesh::invalid_uint));
      66             : 
      67             :   // If we're our own _mesh_sys, we'd better be in charge of
      68             :   // at least one coordinate, and we'd better have the same
      69             :   // FE type for all coordinates we are in charge of
      70             :   libmesh_assert_not_equal_to (mesh_xyz_var, libMesh::invalid_uint);
      71             :   libmesh_assert(!n_x_dofs || context.element_fe_var[_mesh_x_var] ==
      72             :                  context.element_fe_var[mesh_xyz_var]);
      73             :   libmesh_assert(!n_y_dofs || context.element_fe_var[_mesh_y_var] ==
      74             :                  context.element_fe_var[mesh_xyz_var]);
      75             :   libmesh_assert(!n_z_dofs || context.element_fe_var[_mesh_z_var] ==
      76             :                  context.element_fe_var[mesh_xyz_var]);
      77             : 
      78             :   const std::vector<std::vector<Real>>     & psi =
      79             :     context.element_fe_var[mesh_xyz_var]->get_phi();
      80             : 
      81             :   for (auto var : make_range(context.n_vars()))
      82             :     {
      83             :       // Mesh motion only affects time-evolving variables
      84             :       if (this->is_time_evolving(var))
      85             :         continue;
      86             : 
      87             :       // The mesh coordinate variables themselves are Lagrangian,
      88             :       // not Eulerian, and no convective term is desired.
      89             :       if (/*_mesh_sys == this && */
      90             :           (var == _mesh_x_var ||
      91             :            var == _mesh_y_var ||
      92             :            var == _mesh_z_var))
      93             :         continue;
      94             : 
      95             :       // Some of this code currently relies on the assumption that
      96             :       // we can pull mesh coordinate data from our own system
      97             :       if (_mesh_sys != this)
      98             :         libmesh_not_implemented();
      99             : 
     100             :       // This residual should only be called by unsteady solvers:
     101             :       // if the mesh is steady, there's no mesh convection term!
     102             :       UnsteadySolver * unsteady;
     103             :       if (this->time_solver->is_steady())
     104             :         return request_jacobian;
     105             :       else
     106             :         unsteady = cast_ptr<UnsteadySolver*>(this->time_solver.get());
     107             : 
     108             :       const std::vector<Real> & JxW =
     109             :         context.element_fe_var[var]->get_JxW();
     110             : 
     111             :       const std::vector<std::vector<Real>> & phi =
     112             :         context.element_fe_var[var]->get_phi();
     113             : 
     114             :       const std::vector<std::vector<RealGradient>> & dphi =
     115             :         context.element_fe_var[var]->get_dphi();
     116             : 
     117             :       const unsigned int n_u_dofs = context.dof_indices_var[var].size();
     118             : 
     119             :       DenseSubVector<Number> & Fu = *context.elem_subresiduals[var];
     120             :       DenseSubMatrix<Number> & Kuu = *context.elem_subjacobians[var][var];
     121             : 
     122             :       DenseSubMatrix<Number> * Kux = n_x_dofs ?
     123             :         context.elem_subjacobians[var][_mesh_x_var] : nullptr;
     124             :       DenseSubMatrix<Number> * Kuy = n_y_dofs ?
     125             :         context.elem_subjacobians[var][_mesh_y_var] : nullptr;
     126             :       DenseSubMatrix<Number> * Kuz = n_z_dofs ?
     127             :         context.elem_subjacobians[var][_mesh_z_var] : nullptr;
     128             : 
     129             :       std::vector<Real> delta_x(n_x_dofs, 0.);
     130             :       std::vector<Real> delta_y(n_y_dofs, 0.);
     131             :       std::vector<Real> delta_z(n_z_dofs, 0.);
     132             : 
     133             :       for (unsigned int i = 0; i != n_x_dofs; ++i)
     134             :         {
     135             :           unsigned int j = context.dof_indices_var[_mesh_x_var][i];
     136             :           delta_x[i] = libmesh_real(this->current_solution(j)) -
     137             :             libmesh_real(unsteady->old_nonlinear_solution(j));
     138             :         }
     139             : 
     140             :       for (unsigned int i = 0; i != n_y_dofs; ++i)
     141             :         {
     142             :           unsigned int j = context.dof_indices_var[_mesh_y_var][i];
     143             :           delta_y[i] = libmesh_real(this->current_solution(j)) -
     144             :             libmesh_real(unsteady->old_nonlinear_solution(j));
     145             :         }
     146             : 
     147             :       for (unsigned int i = 0; i != n_z_dofs; ++i)
     148             :         {
     149             :           unsigned int j = context.dof_indices_var[_mesh_z_var][i];
     150             :           delta_z[i] = libmesh_real(this->current_solution(j)) -
     151             :             libmesh_real(unsteady->old_nonlinear_solution(j));
     152             :         }
     153             : 
     154             :       for (unsigned int qp = 0; qp != n_qpoints; ++qp)
     155             :         {
     156             :           Gradient grad_u = context.interior_gradient(var, qp);
     157             :           RealGradient convection(0.);
     158             : 
     159             :           for (unsigned int i = 0; i != n_x_dofs; ++i)
     160             :             convection(0) += delta_x[i] * psi[i][qp];
     161             :           for (unsigned int i = 0; i != n_y_dofs; ++i)
     162             :             convection(1) += delta_y[i] * psi[i][qp];
     163             :           for (unsigned int i = 0; i != n_z_dofs; ++i)
     164             :             convection(2) += delta_z[i] * psi[i][qp];
     165             : 
     166             :           for (unsigned int i = 0; i != n_u_dofs; ++i)
     167             :             {
     168             :               Number JxWxPhiI = JxW[qp] * phi[i][qp];
     169             :               Fu(i) += (convection * grad_u) * JxWxPhiI;
     170             :               if (request_jacobian)
     171             :                 {
     172             :                   Number JxWxPhiI = JxW[qp] * phi[i][qp];
     173             :                   for (unsigned int j = 0; j != n_u_dofs; ++j)
     174             :                     Kuu(i,j) += JxWxPhiI * (convection * dphi[j][qp]);
     175             : 
     176             :                   Number JxWxPhiIoverDT = JxWxPhiI/this->deltat;
     177             : 
     178             :                   Number JxWxPhiIxDUDXoverDT = JxWxPhiIoverDT * grad_u(0);
     179             :                   for (unsigned int j = 0; j != n_x_dofs; ++j)
     180             :                     (*Kux)(i,j) += JxWxPhiIxDUDXoverDT * psi[j][qp];
     181             : 
     182             :                   Number JxWxPhiIxDUDYoverDT = JxWxPhiIoverDT * grad_u(1);
     183             :                   for (unsigned int j = 0; j != n_y_dofs; ++j)
     184             :                     (*Kuy)(i,j) += JxWxPhiIxDUDYoverDT * psi[j][qp];
     185             : 
     186             :                   Number JxWxPhiIxDUDZoverDT = JxWxPhiIoverDT * grad_u(2);
     187             :                   for (unsigned int j = 0; j != n_z_dofs; ++j)
     188             :                     (*Kuz)(i,j) += JxWxPhiIxDUDZoverDT * psi[j][qp];
     189             :                 }
     190             :             }
     191             :         }
     192             :     }
     193             : #endif // 0
     194             : 
     195             :   return request_jacobian;
     196             : }
     197             : 
     198             : 
     199             : 
     200    16452128 : bool FEMPhysics::mass_residual (bool request_jacobian,
     201             :                                 DiffContext & c)
     202             : {
     203     1495648 :   FEMContext & context = cast_ref<FEMContext &>(c);
     204             : 
     205     1495648 :   unsigned int n_qpoints = context.get_element_qrule().n_points();
     206             : 
     207    32904256 :   for (auto var : make_range(context.n_vars()))
     208             :     {
     209    17947776 :       if (!this->is_time_evolving(var))
     210           0 :         continue;
     211             : 
     212    14956480 :       FEBase * elem_fe = nullptr;
     213    14956480 :       context.get_element_fe( var, elem_fe );
     214             : 
     215    16452128 :       const std::vector<Real> & JxW = elem_fe->get_JxW();
     216             : 
     217     1495648 :       const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
     218             : 
     219             :       const unsigned int n_dofs = cast_int<unsigned int>
     220     2991296 :         (context.get_dof_indices(var).size());
     221             : 
     222     1495648 :       DenseSubVector<Number> & Fu = context.get_elem_residual(var);
     223     1495648 :       DenseSubMatrix<Number> & Kuu = context.get_elem_jacobian( var, var );
     224             : 
     225    89504800 :       for (unsigned int qp = 0; qp != n_qpoints; ++qp)
     226             :         {
     227           0 :           Number uprime;
     228    73052672 :           context.interior_rate(var, qp, uprime);
     229    79693824 :           const Number JxWxU = JxW[qp] * uprime;
     230   365263360 :           for (unsigned int i = 0; i != n_dofs; ++i)
     231             :             {
     232   345339904 :               Fu(i) -= JxWxU * phi[i][qp];
     233   292210688 :               if (request_jacobian && context.elem_solution_rate_derivative)
     234             :                 {
     235    87408640 :                   const Number JxWxPhiIxDeriv = JxW[qp] * phi[i][qp] *
     236    87408640 :                     context.elem_solution_rate_derivative;
     237    87408640 :                   Kuu(i,i) -= JxWxPhiIxDeriv * phi[i][qp];
     238   218521600 :                   for (unsigned int j = i+1; j < n_dofs; ++j)
     239             :                     {
     240   154951680 :                       const Number Kij = JxWxPhiIxDeriv * phi[j][qp];
     241   131112960 :                       Kuu(i,j) -= Kij;
     242   131112960 :                       Kuu(j,i) -= Kij;
     243             :                     }
     244             :                 }
     245             :             }
     246             :         }
     247             :     }
     248             : 
     249    16452128 :   return request_jacobian;
     250             : }
     251             : 
     252             : } // namespace libMesh